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1 . INTRODUCTION 

SEASAT-A was launched on 26 June, 1978 carrying a complement 
of active and passive microwave instruments designed to sense 
the vector wind stress at the sea surface, the significant wave 
height and other atmospheric and oceanographic surface properties. 
The satellite stopped transmitting after acquiring global data 
from each instrument over a 99 day period. A preliminary over- 
view of the SEASAT-A mission has been reported by Born et al. 
(1979a) . 

The evaluation of the spacecraft and sensor performance must 
still be regarded to be in a preliminary stage. Nevertheless, 
on the basis of intensive analysis of limited geophysical series, 
derived from sensor records by a set of geophysical retrieval 
algorithms developed before launch, it is apparent that the design 
goals for measurement accuracy will be met, at least for those 
sensors of interest here (Born et al., 1979b). The research pro- 
gram, of which this study is a part, is concerned primarily with 
radar backscatter measurements from the scatterometer (SASS) 
and wind vectors (affected by aliases) derived therefrom; with 
scalar winds derived from the passive microwave measurements by 
the multichannel radiometer (SMMR) ; and with significant wave 
height measurements from the radar altimeter (ALT) . 

Quantitative wind and surface wave data, available on a 
global basis, shall have great potential for improving our 
knowledge of and ability to predict weather and sea state through 
the improved specification of initial conditions for numerical 
weather and ocean models. However, because of the novel nature 
of such synoptic data it is not known exactly how the data should 
be used to achieve maximum impact, how much improvement in 
environmental forecasts can be expected, or how best to exploit 
SEASAT measurement capabilities in a future operational system. 

Following the refinement and evaluation of algorithms for 
retrieval of geophysical data from seasat sensor data, global 
wind and wave data sets will be made available to the scientific 
community for studies directed toward resolution of the above 
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questions. Prior to the availability of real' data sets, however, 
answers to the above questions can be sought through simulation 
experiments carried out with general circulation (GCM) models 
and coupled ocean-GCM models* Such experiments should also lead 
to the development of optimum techniques for the assimilation of 
real SEAS AT -A data into numerical weather and sea state fore- 
casting models. 

The very first SEASAT simulation experiments began in 1975 
and were conducted at the Goddard Institute of Space Studies 
(GIoS) as a collaborative effort involving several groups, as 
reported by Cardone (1977) and Haiem et al. (1976). This study 
describes more recent results of the collaborative program with 
particular emphasis on three areas in which Oceanweathet Inc. 
contributed significantly: (1) the design and analysis of SEASAT 
simulation studies in which the error structure of conventional 
analyses and’ forecasts is modeled realistically; (2) the develop- 
ment and implementation on the GMSF computer of a global spectral 
ocean wave model (GSOWM) which may be run in tandem with the 
GMSF GCM to assimilate SEASAT data and assess forecast sea state- 
impacts; (3) the design of algorithms for the assimilation of 
SASS wind data into the GMSF GCM and for the utilization of real • 
SASS wind data and ALT wave height data in a coupled GCM-GSOWM. 
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2, SEASAT SIMULATION STUDIES 
2 . 1 Background 

A general experimental design for SEASAT simulation studies 
has been formulated through the efforts of a group of scientists 
engaged in general circulation modelling, numerical weather 
prediction, and SEASAT-related research. The group first met at 
the National Center for Atmospheric Research (NCAR) in February 
of 1975# and produced an overall plan of action that suggested 
near-term, intermediate, and long-range goals. The group now 
serves as an ad hoc SEASAT simulation studies steering group and 
meets periodically to review results of the simulation studies 
and to participate in the refinement and evaluation of the experi- 
ment plan. 

The first experiments suggested were of the so-called iden- 
tical-twin type, in which a specific general circulation model 
provides both the "real" atmosphere reference data and the 
simulated satellite data, and is used to assimilate the satellite 
data to study the impact of such assimilation on subsequent fore- 
casts. Those experiments were to simulate SEASAT-A, and a pos- 
sible multiple SEASAT operational system, and were to involve 
straightforward insertion of simulated SEASAT wind measurements 
both alone and with simulated temperature-sounding data, 

Experiments designed largely along those lines were per- 
formed at GISS with the GISS GCM during 1976. The results of 
those experiments have been reported in detail by Halem et al. 

(1976) and Cardone (1977) and McCandless and Cardone (1976). The 
results of those first experiments were mildly encouraging but 
were regarded with great caution because of the well known short- 
comings of the experimental methods employed. In particular, the 
experiment was biased toward optimistic results because the iden- 
tical 3CM had been used in all three runs that comprised the 
simulation study. Also, the random nature of the error specifi- 
cation used to model operational (not enhanced by SEASAT) analyses 
did not realistically simulate the error structure in real analyses 
based upon the conventional observation network. Finally, the 
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asynoptic nature of SEASAT data was not modelled and the method 
used to assimilate the data was very crude. 

Following those experiments, a detailed plan for SEASAT 
simulation experiments was drafted in collaboration with key 
members of the steering group. The plan called for a major re- 
design of the simulation impact tests from the methodology that 
had been employed in most prior satellite simulation studies. 

In particular it recommended the simulation in great detail, of 
the conventional meteorological observing network and the conven- 
tional methods for objectively analysing the observations from 
such network to produce initial states for numerical weather fore- 
casting, with the goal of realistically simulating the error 
characteristics of conventional meteorological data fields and, 
thence, the forecast error growth rates found in nature. The 
realistic control states so produced were to be used to assimilate 
simulated SEASAT data asymptotically. 

The plan called for relatively simple SEASAT simulation ex- 
periments at first, e.g. the direct insertion of perfect SEASAT- 
derived winds, following a seasat-A orbit, at the lowest active 
level of the GISS GCM, i.e. Level 9. Experiments to follow were 
to include simulation of the improvements in marine surface 
pressure to be expected to accompany the SEASAT winds; insertion 
of all satellite data (e.g. VTPR, GOBS winds); and simulation of 
possible future configurations of satellite orbits and sensors. 
Also, it was proposed that a sea-state forecast model be coupled 
to the GCM, to allow assessment of the impact on wave forecasts 
to follow from indicated levels of impact on weather forecasts. 

The evolution of simulation experiments was to proceed to the use 
of more sophisticated assimilation methods, and ultimately to 
experiments with real SEASAT-A data. 

2.2 Summary of Results 

Significant progress has been made to date toward the 
achievement of the objectives stated above. A set of simulation 
studies has been completed, using the GMSF GCM, which modelled 
the error structure in conventional analyses and forecasts far 
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more realistically than has been done in prior studies. It has 
been possible to demonstrate, rather convincingly, that against 
such control analyses and forecasts, SEASAT surface wind data 
have the potential to improve numerical weather forecasts signif- 
icantly over extra tropical ocean and land areas. 

The detailed results of the realistic simulations have been 
documented and submitted for publication in revised form as a 
paper entitled: 

"Observing Systems Simulation and Potential impact 
of Marine Surface Wind Data on Numerical Weather 
Prediction" by M, Cane, V, Cardone , M. Halem, 

I. Halberstam, J. Ulrich (Submitted to Mon, Wea. 

Rev. ) 

A copy of the paper is attached to this report as Appendix C. 

Only a summary of results is given here. 

A simulation experiment requires four elements: [13 a "nature 
run", i.e. a month- long integration of the GCM from an arbitrary 
initial state, to provide a complete record of the state of the 
model atmosphere, which is used noth to fabricate "observational" 
reports and to evaluate analyses and forecasts; [2] a control 
assimilation that is like an operational forecast-analysis cycle 
based on conventional observations, except that the data used to 
produce the analyses and fields are in fact fabricated from the 
nature run; [33 a SEASAT assimilation, differing from the control 
assimilation in incorporating "SEASAT data" (likewise fabricated 
from the nature run) in the forecast-analysis cycle* [4] forecasts 
produced from both initial conditions [2] 6 [3] , which are compared 
to nature and inter se to provide an assessment of the anticipated 
impact of the SEASAT data. 

In the most recent series of simulations, a 30-day history 
run with the GMSF GCM was used to fabricate simulated observations 
at the times and locations of the individual conventional reports 
(surface, radiosonde, and ship) actually received during February 
1976. Those fabricated observations, suitably degraded for instru- 
ment and sampling errors, were then used to create analysed fields 
on the GCM grid in an analysis-forecast cycle like those in use 
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at major meteorological centers. The control fields so produced 
are much more representative of actual analyses than those pro- 
duced by perturbing initial states with random errors. The 
forecast error growth in five simulated 72-hour forecasts from 
the control states was gr&tifyingly similar to that found in 
operational numerical forecasts. 

A series of experiments has been conducted, each simulating 
the addition of surface winds derived from SEASAT-A to the control 
run. Xn all those experiments, the SBASAT winds were fabricated 
directly for the GCM grid points intercepted by the SEASAT-A 
scatterometer swath, and representative of the lowest level of 
the GCM (about 950 mb) . Two of the experiments assumed error- 
free SEASAT winds: the first assimilated wind data asynoptically 
by the direct insertion method (PW-DIM) ; the second employed the 
successive correction method (PW-SCM) . An experiment was also 
performed which simulated the addition of perfect sounder derived 
temperatures by the SCM (PT-JCM) to the control, so that the 
relative impact of surface wind data compared to sounder data 
(both error free) could be assessed. The results of forecast 
experiments from each of these experiments are compared briefly 
here. (A more detailed discussion is contained in Appendix C) . 

In each of the experiments, 72-h forecasts were made from 
OOQOGMT on 5 Feb., 10 Feb., 15 Feb., 20 Feb., and 25 Feb., and 
compared to control forecasts and the verifying states from the 
nature run. Comparisons were made in terms of growth of rms 
error in level 9 zonal wind speed (U 9 ) , surface pressure (P g ) , 
and level 5 (around 500 mb) zonal wind speed (Ug) . In addition, 
subjective evaluation of isobaric sea level and upper level 
contour maps was performed for several cases. 

Table 1 displays rms errors relative to nature in P g and 
U 5 for individual forecasts from the control, PW-DIM, PW-SCM and 
PT-SCM initial states, with the statistics sorted by day and re- 
gion. Statistics for Ug and for other regions (Eurasia, North 
Atlantic, South America) were similarly prepared and studied. 

The errors for day 0 represent the errors in the initial states. 
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Because of the small sample size, there was no attempt to establish 
the statistical significance of the average differences between 
the control and satellite forecast? however, as an example, the 
average growth of error is shown graphically in Pig. 1 for P g 
so that relative differences in forecast errors can be assessed 
subjectively. 

The results of the PW-DXM experiment generally confirmed the 
expectation that direct insertion of level 9 winds would have no 
significant impact on analyses and forecasts. The only significant 
response seen is about a 10# average Improvement in the level 9 
zonal (and meridional) wind component in the initial states over 
oceanic regions. However, errors in the initial sea level pres- 
sures suggest a tendency for the PW-DXM surface pressure field to 
be degraded slightly by the wind data assimilation, especially in 
the tropics. The isobaric analysis examined in the tropics, where 
mean gradients are weak, revealed that small scale variations had 
been added to the PW-DXM fields. The variations are probably 
associated with shock attending the DIM assimilation. 

A comparison of forecast rms errors between Control and PW-DXM 
forecasts shows two interesting effects. Primarily, overall mean 
impacts are small in all variables and regions. The improvement 
in Ug fields is usually lost with the first forecast day. How- 
ever, subjective examination of sea level pressure charts did 
indicate a few noticeable impacts in individual forecasts, usually 
shown as improved placement and intensity of extratropical cyclones 
over oceanic areas. 

The assimilation of error free level 9 winds by successive 
correction (PW-SCM) resulted in much larger and more consistent 
positive impacts on analyses and forecasts. Reductions of rms 
Ug errors relative to the Control or PW-DIM analyses averaged 
about 30#. Initial surface pressure errors, which were larger 
than the Control for the PW-DIM analyses, are slightly lower than 
Control errors. The improvements in the Ug and P g analyses are 
consistent from region to region and from day to day (see Fig. 1). 
Average improvements in the Ug PW-SCM fields are 30#, except over 
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the tropical oceans and South America, where only slight impacts 
were found. 

The forecasts made from error free wind SCM analyses are 
generally more skillful than corresponding Control and PW-DIM 
forecasts by margins equivalent to relative differences in the 
initial states. In contrast to the PW-DIM forecasts, which 
displayed occasional positive impacts, the forecasts of Ug and P D 
are improved over Control and DIM forecasts in all cases and all 
regions. On the other land, at level 5, Improvements in the PW-SCM 
initial states did not always lead to improvements in the forecasts. 

A subjective evaluation of the improvements showed that the 
predominant effect of the assimilation of SASS data was to greatly 
reduce the large and spatially coherent errors which characterized 
the control level 9 over ocean windfield. That improvement in 
turn led to substantial improvements in the 3 day forecast level 9 
wind field and surface pressure field both over ocean areas and 
over downstream continental regions. It was apparent also that 
the impacts favored the eastern North Pacific Ocean, where Seasat 
SCM forecast pressure analyses revealed consistently better fore- 
casts of the intensity and placement of extratropical cyclones. 

This feature of the simulations may be related to the poor cov- 
erate of conventional data in the North Pacific, relative to the 
North Atlantic and continental regions, and to the fact that the 
simulated Seasat orbit favored the insertion of simulated Seasat 
wind data over the Pacific during tne 6-h period prior to 0000 
GMT, which was the initial time for all forecast simulations. 

The PT-SCM experiment was conducted to provide an indication 
of the relative impact of surface wind data compared to sounder 
data, when both sets of observations could be considered to be 
error free. Most previous simulation studies have in fact dealt 
with sounders. While the earlier studies have tended to be quite 
optimistic regarding the potential impact of sounder derived 
temperatures on NWP, much of the optimism probably stemmed from 
unrealistic Control experiments and the assumption of sounder 
errors much lower than have been attained operationally. Recent 
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studies with real data, have shown that sounders can have iimall, 
but statistically significant positive impact on NWP. The im- 
pact appears to be highly sensitive to the quantity of data and 
the assimilation method. For Nimbus sounder data, assimilated 
time-continuous ly by SCM, Ghil et al. (1979) find 5 - 1 % positive 
impacts in rms measures of forecast skill over North America. 

The results of the simulated perfect Nimbus souding SCM 
experiment can be seen in Table 2 and Fig. 1* In general the 
impacts found were comparable to those of the PW-SCM experiment. 
The only difference evident is that the sounding data is slightly 
more effective than the surface wind data in controlling growth 
of level 5 forecast wind errors. Within the context of our 
experiment, the two data types appear to have equivalent value 
in reducing surface pressure errors. The perfect sounding pro- 
duced about a 1555 improvement in sea level pressure (rms) fore- 
cast errors. This appears to be a more reasonable upper limit 
to impacts to be expected from satellite data than have been 
predicted in less realistic OSSE's. 

It is reasonable to assume that the idealized error free 
wind experiments overestimate the positive impact to be expected 
from real Seasat-1 data. However, a limited SCM experiment was 
performed to assess the effect of nominal errors on the satellite 
winds (±2m/s in speed, ±20° in direction, normally distributed) . 
The results of that experiment are compared to the error free 
experiment in Fig. 2. Evidently for the fairly dense distribution 
of remotely sensed winds, the SCM assimilation is very effective 
at removing errors that are uncorrelated. Actual scatteroraeter 
wind data is likely to have a more complicated error structure, 
but may significantly impact forecasts if the errors are as small 
overall as the above, and if a good enough assimilation scheme 
is used. 

A comparison of the PW-SCM and FT-SCM experiments suggests 
that surface wind data has the same potential impact as tempera- 
ture sounders, when both sets of observations are error free, es- 
pecially over the downstream of the eastern North Pacific and 
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North Atlantic basins. Cane et al. (1979) contend that indirect 
support for this result is provided by Blackmon, et al. (1979). 

They show from observational data that over the eastern sides of 
the Northern Hemisphere oceans the 500 mb height is much more 
strongly correlated with 1000 mb height than with 1000-500 mb 
thickness and infer that surface (lOOOmb) data would play an 
important role in the determination of mid tropospheric structure. 
Satellite sounding data are of greatest potential over continents, 
where Blackmon et al. show that 500 mb height is more strongly 
correlated with 1000-500 thickness than with 1000 mb height. How- 
ever, over continents, sounding data are largely redundant with 
conventional radiosondes . 

The mechanism which Cane et al, (1979) hypothesize as respon- 
sible for the relatively large impacts of simulated Seasat data 
emerges as follows. Given the essentially barotropic nature of 
the atmosphere over the eastern North Pacific and North Atlantic, 
improvements in the surface (pressure) fields can significantly 
impact tropospheric analyses and forecasts. In terms of geo- 
strophic adjustment theory, the scale of the wind data assimilated 
is small compared to the (barotropic) radius of deformation. There 
fore it is reasonable to expect wind information to be retained 
and the mass field to adjust, while temperature data will tend to 
be radiated away as gravity waves. Further, from a statistical 
point of view, there is more information in a small orbit segment 
of wind data than in sounder data of the same size because the 
correlation scale of winds is smaller. This suggests that wind 
data can benefit more from the dense coverage that a satellite 
provides . In addition , even the preliminary evaluation of actual 
&ASS data suggests that scatterometer wind errors are a smaller 
fraction of marine boundary layer wind analysis error than are 
sounder temperature errors compared to temperature analysis. 

The results of these idealized impact studies, coupled with 
the fact that actual Seasat- A SASS marine wind data have an 
accuracy close to nominal specifications suggest that studies 
involving real SASS global data sets should be undertaken and that 
serious consideration be given to such data in the design of an 
optimum global observing system. 
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3. FUNCTIONAL DESCRIPTION AND DOCUMENTATION OF GMSF GLOBAL 

SPECTRAL OCEAN WAVE MODEL (GSOWM) 

To implement a global spectral ocean wave model (hereafter 
GSOWh/ usable on the Amdahl 470 at Goddard Modelling and 
Simulation Facility (hereafter GMSF) and readily portable to any 
IBM 370 of roughly the same size, three programs were written 
and tested. Programs PRELIM (to print grid coordinates and 
propagate tables, and write propagate tables to disk) and 
ICEDECK (to write land-sea tables to disk) need be run once each, 
at installation time; program PAXXPAXX contains the growth and 
propagation modules, and is run for every time segment it is 
proposed to simulate. 

The auxiliary programs DUMP 1 3 , DUMQ13, DUMP 14 permit 
retrieval of the records written to the output datasets by 
PAXXPAXX. 

3.1 PRELIMS. This program, for computing propagation 
coefficients on a truncated latitude-longitude grid (the polar 
caps are deleted) , derives from a program submitted to FNWC in 
1977 for a full latitude-longitude grid (including the poles) 
with improvements discovered while adapting the program to a 
transverse Mercator grid for NOAA/AOML/SAIL in 1978. The prin-? 
cipal novelty introduced in 1978 was the uniform use of a triangle 
on the earth as a basis for computing the outbound propagate 
table: in 1977 rectangles were used in general, but triangles 
where the longitudinal spacing changed. The principal novelty 
introduced here is the use of a distinct scratch file for each 
propagable frequency (14 in the model submitted : see §3.1.1) which 
remarkably simplifies the housekeeping operations in transposing 
the propagate table ( § 3 . 1 . 5 . 7 ) . 

3.1.1 Data definitions 

13 Date, set containing the transposed propagate table: 

79 tracks. This dataset is referenced every time 
a wave field is simulated; if disk storage is labile, 
a backup tape is advisable. 
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15 through 28. Datasets containing the inbound propagate 

tables for the 14 propagable frequencies: 6 tracks 
each. After successful execution of PRELIMS, these 
can be decatalogued. 

3.1.2 Grid numbering system . 

Grid points are numbered by the subscripts (LONG, LAT) , 

1_<LONG<76 , 1 <LAT<73 . LONG increases from west to east; LAT 
increases from south to north. Grid points may also be ref- 
erenced by the single subscript KPOINT w LQNG+76* (LAT-1 ) . 

Denote true latitude (degrees, north positive) by ALAT; true 
longitude (degrees, east positive) by ALONG. Then 
ALAT = -74.0+4. Of LAT, 

ALONG - AMOD (80.0+2. 5+MOD (LAT , 2) +5 . O+LONG, 360.0)-180.0. 

On the globe, points with ALONG = 73, 74, 75, 76 over i.e. points 
with ALONG ■ 1, 2, 3, 4. This unusual choice of prime meridian 
obviates a coding exception to prevent waves from propagating 
across the Isthmus of Panama: the land-sea table handles the 
exception. 

3.1.3 Input parameters . 

In the program as submitted there are none. In a generalized 
program accommodating arbitrary numbers of frequency and direction 
bands, the numbers 14 and 20 occurring in DIMENSION statements 
would be PARAMETER variables; the variables FREQ1 and FREQ2 
would be input parameters. Then compute WEDGE as 186 divided 
by the number of directions; compute RAT = (FREQ2/FREQ1 ) ** ( 1 . / (NF-1 . 5) ) , 
where NF is the number of frequencies; compute DELT as discussed 
in § 3. 1 . 5. 4. 

3.1.4 Common blocks . 

N.B. All variables in blank common are local and volatile. 

COMMON 

$ /YPARAM/ LLAT , LLONG, LPNT1, LPNT2, FREQC14 ), FREQl, FREQ2, SPACE(2 ) 

$ , WEDGE, DIREC(20),DELT 

$ / YSMORT / TLAT(3,73),TLONG(2,76) 

$ /YYOUTB/ JTABLE,TABLEC4, 20,73) 

$/YYTRIG/ DLAT, COSS, DANGLE, CLAT, TRAVEL, RAD, CTR, STR, 

$ OMIN, 017, RAT (3/ 73) 
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YPARAM: 

LLAT number of parallels, i.e. 73 

LLONG number of meridians, i.e. 76 

LPNT1 number of grid points, without deduction for land, but 
with duplicates excluded: 73x72 « 5256 
LPNT2 gross number of grid points; 73x76 * 5548 
FREQ nominal frequencies of the 14 bands 
FREQ1 nominal frequency of lowest band, i.e. .04 hertz 
FREQ2 nominal frequency of highest band, i.e. .20 hertz 
SPACE approximate maximum and minimum values of grid spacing, nm 
WEDGE half of the angular bandwidth, i.e. 9° 

DIREC array of 20 directional bands, in degrees, clockwise 

from north: DIREC(I) = 18*1-9. Note : in 1977 we found 

that propagation along a meridian or the equator entailed 
a vexatious exception to the general code; directional 
bands are now chosen to miss the cardinal points. 

DELT time step used in propagation, i.e. 3 hours. See §3. 1.5. 3. 

YSHORT: 

TLAT functions of grid latitude. See §3. 1.5. 2. 

TLONG functions of grid longitude. See S3. 1.5.2. 

YYOUTB : 

JTABLE grid point numbers and angular band numbers in outbound 
propagation table. See §3.1 .5.5. 

TABLE energy fractions in outbound propagation table. See §3. 1.5. 5. 
YYTRIG: 

DLAT colatitudes, in radians 
COSS cosines of the 20 directions 

DANGLE the 20 directions, in radians, reduced to the range 
(-it, -Hr) 

CLAT cosines and sines of the colatitudes 

TRAVEL distance covered by waves of each frequency, at their 
group velocity, in time DELT, in nm. See §3. 1.5.4. 

RAD 1 80/it 

CTR cosines of TRAVEL 


I / 

STR sines of TRAVEL 

OMIN tt/10800 

017 tr/1 80 

RAT correction factors for convergence of meridians. See S3. 1.5. 4. 

3.1.5 Subprograms . 

3. 1.5.1 MAIN. Defines common blocks, calls GGRID, BANDS, and 
TRIG (once each) , calls OUTBND followed by INBND (once 
for each frequency) , calls TRANSP once, and writes a 
file mark on dataset 13. (This step was inserted when 
coding for SAIL, where all permanent files are tapes) . 

3. 1.5. 2 GGRID. The principal products of the program are 
arrays TLAT (referenced by subroutine TRIG) , and JLAT, 
JLONG, ALAT, ALONG (used to print table of coordi- 
nates) . JLAT, JLONG, ALAT, ALONG become obvious upon 
inspection of the printed output; TLAT, TLONG require 
explanation. TLONG (1,1) is the true longitude of 

the grid points (1,1), (1,3), etc.; TLONG (2,1), of 
the points (1,2), (1,4), etc. To understand TLAT, 
draw the 73 parallels, as well as zigzag lines con- 
necting grid points on adjacent parallels, yielding 
a mesh of isosceles triangles: the point (ALAT, ALONG) 
has the 6 neighbors (ALAT, ALONG+5) , (ALAT±2 , ALONG+2.5). 
TLAT (1,1) true latitude corresponding to grid latitude 1. 
TLAT (2,1) base of isosceles triangles, i.e. 5° of 
longitude, in nm. 

TLAT (3,1) approximate mean value (flat earth computa- 
tion) of legs of isosceles triangles. 

3. 1.5.3 BANDS. Computes the midpoints of the 20 angular bands 
and the nominal frequencies of the 14 frequency bands. 

The choice of frequency-direction bands (or bins) for 
a SOWM has been more arbitrary than scientific. Uji 

6 Isozaki have used bands uniformly spaced in period. 
Pierson in principle uses elementary bands of width 
1/180 hertz (to agree with an ancient Tukey spectral 
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analysis of wave records at British weather ships) ; 
at higher frequencies, he groups 2, 3, or 4 elementary 
hands into one working band. If wave spectra are 
approximately similar, then spectral resolution is 
equalized at all sea states by spacing the frequencies 
in geometric progression. A little trial showed that 
a ratio of 5 between the highest and lowest bands was 
adequate, but that the absolute frequency of the 
lowest band may need to be shifted according to the 
wave regime under study. The desirable number of 
angular bands is also uncertain. Gelci in his early 
experiments used 8 bins. Pierson et al. prepared a 
1 2-angular~band model for FNWC, and concluded from 
their numerical tests that 24 bands would yield much 
better representation of swells. We used 24 bands 
in the limited-area SOWM prepared for SAIL. If the 
purpose of banding is a faithful modelling of prop- 
agation, then the termini of group-velocity vectors 
corresponding to a frequency-direction band should 
approximate the vertices of a square, leading to an 
angular bandwidth (in radians) of log e • 

With the value of RAT used in BANDS, the angular band- 
width would be 7.38°: in round numbers 7.5°, or 48 
bands. We are aware of no numerical experiments with 
angular resolution this fine, perhaps because very 
few series of directional spectra have been acquired 
to calibrate such a model. 

3. 1.5. 4 TRIG. Computes the great-circle distance traveled in 
one time step at each nominal frequency, and trigono- 
metric functions of these arcs, colatitudes, and 
directions of propagation. These calculations are in 
double precision because we will be solving an ill- 
conditional spherical triangle: two long sides (co- 
latitudes of origin and terminus of propagation vector) 
and one short side (distance traveled, which is less 
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than 2° . ) The time step DELT is chosen as the largest 
multiple of 1 hour in which the .04 hertz waves 
travel less than 120 nm. TRIG also computes the area 
ratio between adjacent parallels! each grid point 
is treated as a box 120 nm x (300 cos ALAT) nm# and 
the quantity conserved in propagation is the wave 
variance multiplied by the area of the box. 

3. 1.5. 5 OUTBND. Executes steps 1, 2, 3 of S3. 1.6. The number- 
ing of the eight quadrant indicators is: 
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Triangles III and IV occur only at latitudes greater 
than ±66° and frequencies of .04 and .045496 hertz. 

In numbering the vertices of the triangles, the upper 
sign affixed to ALAT applies to the NE and NW quadrants 
the upper sign affixed to ALONG applies to the NE and 
SE quadrants. Waves start at the point (ALAT, ALONG) . 


Triangle I , vertex 1 : 

2 : 

3: 

Triangle II, vertex 1 : 

2 : 

3: 


ALAT ±2.0, ALONG ±2.5 
ALAT ±2.0, ALONG *2.5 
ALAT, ALONG 
ALAT, ALONG ±5.0 
ALAT ±2.0, ALONG ±2.5 
ALAT, ALONG 
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3. 1.5. 6 


3. 1.5.7 


Triangle III, vertex 1: 

2; 

3: 

Triangle IV, vertex 1 s 

2 ; 

3: 


ALAT ±2.0, ALONG ±7,5 
ALAT, ALONG ±5.0 
ALAT ±2.0, ALONG ±2.5 
ALAT, ALONG ±10,0 
ALAT ±2.0, ALONG ±7.5 
ALAT, ALONG ±5.0 


INBND. Executes step 4 of 53,1.6. 


TRANSP. Executes step 5 of §3.1.6, In order to econ- 
omize storage by using the type INTEGER* 2 , the present 
code differs from previous transpose routines in that 
the arrays JTAB4 and TAB 4 are not equivalenced. 
Contents of JTAB4 and TAB 4 : 


IFR indexes frequency 
IDIR indexes direction 
LAT is grid latitude 

JTAB4 (1 , IPR, IDIR, LAT) is an integer in the range 1 
to 6 (on the average 4) : the number of parts in the 
propagation formula. 

JTAB4 (2*J, IFR, IDIR, LAT) is the transmitting grid 
point number minus the receiving grid point number. 
Used to enter the land-sea table and verify that 
the transmitting grid point is water before propa- 
gating. 

JTAP4 (2*J+1 , IFR, IDIR, LAT) is the relative address 
of the transmitting band from the receiving band, 
biased by -21280. 

TAB4 (J,IFR, IDIR, LAT) is the fraction of energy in 
the transmitting band to be moved to the receiving 
band. 


3.1.6 Construction of propagate table . 

1 . Start a unit component (one direction-frequency band) , multi- 
plied by the area of a rectangle centered on a grid point, 
and propagate it from that point at its proper group velocity 
and along an arc of a great circle for a time step. In 
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general the receiving point will be in the interior of a 
triangle bounded by grid points, and the true azimuth at 
the receiving point will not be one of the 20 azimuths used 
in CMPE27. 

2. Propagate a suitable fraction of energy into each of the 3 
grid points, corners of that triangle. Turn a fraction of 
one of those fractions through ±15° (toward the equator) to 
conserve the mean azimuth at the receiving point. The search 
for a fraction to turn is made over vertices 1, 2, 3 in order, 

3. Divide each of the 4 energies found in (2) , which sum 
exactly to the energy transmitted in (1), by the area of the 
rectangle centered on the receiving grid point, thus recon- 
verting them to spectral components. 

4. The result of step (3) is an outbound propagation matrix, 
which tells where the energy at time step N goes at step 
N+1. Xt is now necessary to transpose this matrix, yielding 
an inbound propagation matrix containing the multipliers to 
be applied to the several transmitting grid points and 
directional bands at time step N-1 to o'rtain one band at one 
receiving grid point at step N. The number of multipliers 
averages 4: it varies from 1 to 6, in the model submitted. 

The numbers 0 and 7 have occurred with other grids. 

5. The result of step (4) is an inbound propagation matrix 
for each latitude and each propagable frequency: a further 
transpos ition brings together the coefficients for each 
latitude, yielding the propagate tables JTAB4 and TAB4 . 

3.2 XCEDECK . Turns the data statement (written for conve- 
nience in punching) into a printed land-sea table and an array 

stored in DD11. The convention is 1 - land, 2 = sea. 

3.3 PAXXPAXX. 

3.3*1 Data definitions 

10 Dataset containing the wind field, at 4-hour intervals: 

3k tracks per field. The length of a logical record 
is 11096 words, dimensioned (2,76,73). 
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Word (1 ,LONG,LAT) contains wind speed , m/sec, at the 
grid point (LONG,LAT), 

Word ( 2 , IjONG , LAT) contains wind direction, to which 
the wind is blowing, in degrees, clockwise from south, 

11 Dataset containing the land-sea table written by 
1CEDEGK. 1 track. 

12 Dataset containing the propagate table, i.e. DD13 
written by TRANSP (83.1,5.7) 79 tracks. 

13 Principal output file (tape)*, 3.13 megabytes for each 
time step declared by KSTEP3. For format of records 
see §3.3.3. 

14 Supplementary output file (disk); 2 tracks per time 
step. For format of records see §3.3.3. 

17 Scratch file (disk); 241 tracks. 

18 Scratch file (disk): 241 tracks. 

3.3.2 input parameters (in namelist $WI1AT) . 

NAMELIST /WHAT/ LWIND, LLAND, LPROP, LOUT , LARCH1, LARCH2, 

% LRST1, LRST2, LSCR1, L5GR2, KSTEPl, KSTEP2, KSTEP3, YMDH, LSKIP, KBT, 

$ CA, CB, RUN I D 

N.B. The parameters LWIND, LLAND, LPROP, LOOT, LARCH 1 , LARCH 2 , 
LRST1, LRST2 , LSCR1 , LSCR2, were included in case the default 
values of DDs should be unavailable at an installation. Do not 
use them. 

KSTEPl Number of first time step. For dead start, KSTEP1 = 1; 
to continue a previous simulation, KSTEPl (new run) « 
KSTEP2 (previous run) + 1 . 

KSTEP2 Number of last time step. The length of a simulation is 
(KSTEP2-KSTEP1+1) steps, i.e. 3* (KSTO2-KSTEP1+1) hours. 
KSTEP3 Frequency of full output. The tape LARCH 1 is written 

when MOD (TSTEP ,KSTEP3) = 0. Plausible values for KSTEP3 
are 2,4,8; these archive wave spectra every 6,12,24 hours. 
YMDH Date and time corresponding to ISTEP =» 0 , supplied as an 
8-digit integer YYMMDDHH: for example, 79092511. 

LSKIP Number of logical records on DD10 to be skipped before 
reading the first wind field wanted, (default = 0) . 
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KBT The default value, KBT * 1, produces a hump (as used by 
Inoue at wU^cosO/g * .031? to remove the hump, set 
KBT * 0. 

CA The multiplicative constant in the linear growth mechanism. 
The default value is 1.36E-9; any alternative value (for 
sensitivity studies) must be in the deck $WHAT. 

GB The multiplicative constant in the exponential growth 
mechanism. The default value is ,1066; the value ,0231 
is believed to be physically correct. 

RUNID A 9-digit integer serving as a mnemonic to Identify a 
simulation run. 

3.3.3 Format of output data sets . 

DD 1 3 : each time step writes 147 logical records followed by a 

file mark. Logical record 1 contains 3 words, INTEGER* 4 , 
of identification? the time step, the current date and 
time, and the run ID. Logical record (2*LAT) contains 
76 words, REAL *4 , the largest energy (wave variance, 
m ) in one band at the point (LONG,LAT) . Words 1, 2, 

75, 76 are meaningless. Logical record (2+LAT+1) con- 
tains 21280 words, INTEGER* 2, representing an array 

dimensioned (14,20,76). Word (IFREQ, IDIR,LONG) contains 
z 

the variance (m ) in the spectral band (IFREQ, 1DIR) at 

the point (LONG,LAT) , divided by the largest variance 

l s 

in any band at that point, scaled by 2 -1 , rounded and 

converted to integer. 

DD14: each time step writes one logical record, containing 4 

words, XNTEGER+4, of ide u ,ification (array JFRQNT) 
followed by 5548 words, REAL* 4 , dimensioned (76,73), 
containing the significant height, in meters, at the 
point ( LONG, L AT) . 

N.B.: when JFRONT(4) = 0, a full spectral record is 
written to DD13. 

DD18: same format as DD13, but no file mark. Before continuing 

a simulation, copy DD18 to a backup (disk or tape). 
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3.3.4 Common blocks 

COMMON 

$ /XXXXIO/ LWIND, LLAND, LPROP, LOUT, LARCHl, LARCH2, LRSTl, LRST2, 

$ LSCR1, LSCR2, RUNID 

$/XXTIME/ KSTEP1,KSTEP2,KSTEP3, I STEP, DELT, DELTAH, YMDH, ZMDH, 

$ LSKIP, KBT 

$ / XSPEC1/ SPOLDC280, 76, 3) 

$/XSPEC2/ SPNEWC280,76) 

$/XPARAM/ CA, CB,FREQ1,FREQ2,WEDGE,DIREC(20),COSS(2,10), 

$ DELOG, OMEGAC2,l4),FREQCl4),DOMCl4),RBW(14),OMM4Cl4), 

$ 0M3DDC14),0M116C1 y,OM25Cl4),DOMM4Cl4),OMKCl4) 

$/XXGRlD/ LA T , LONG, KPOINT 
$/XXLAND/ LANSEAC76, 73) 

$/XXWIND/ WINDC2, 76, 73) 

$ / XXXH13/ FRONTC4),H13C76,73) 

$/XXPROP/ JPROPC13, 280,73)/PROP(6,280,73) 

$/XXPACK/ EEE(76),SPEC16C280,76) 

3. 3.4.1 XXXXIO 
LWIND = 1 0 
LLAND = 1 1 
LPROP = 1 2 

LOUT * 6 

LARCHl = 13 
LARCH 2 = 14 

LRSTl = 15. Not used. 

LRST2 = 16. Not used. 

LSCR1 = 1 7 
LSCR2 = 1 8 
RUNID See §3.3.2. 

3 . 3 . 4 . 2 XXTIME 

KSTEP1 , KSTEP2 , KSTEP3 , YMDH, LSKIP, KBT See §3.3.2. 

DELT: the growth time step {h the propagation time step) in 

seconds, i.e. 5400 seconds. 

DELTAH: the propagation step in hours, i.e. 3 hours. 

ISTEP: a do- loop index, running from KSTEP1 to KSTEP2. 

ZMDH: current value of date and time, in same format as YMDH 

( see §3.3.2). At time ISTEP = 0, ZMDH is set equal to 
YMDH; at each subsequent time step, ZMDH is updated by 
adding DELTAH and reducing as necessary modulo 24 
hours, 28,29,30,31 days, 12 months, and 100 years. 


25 


3. 3. 4. 3 XSPEC1 

SPOLD The spectra before propagation, at the three parallels 
LAT, LAT± 1 . 

3. 3. 4. 4 XSPEC2 

SPNEW The spectra at the parallel currently being computed. 

3 . 3 . 4 . 5 XPARAM 
CA,CB See 03.3.2. 

FREQ 1 ,FREQ2, WEDGE, DIREC See 03.1.4. 

COSS Array of cosines and sines of multiples of 18°. 

DELOG - ALOG (FREQ2/FREQ1 ) /1 2 , 5 , i.e. the ratio between frequencies 
of successive bands. 

OMEGA Array of radian frequencies. OMEGA(1,I) is the lower 
limit of band X; OMEGA (2, I) is the nominal frequency, 

FREQ Array of nominal frequencies, in hertz. 

DOM Array of bandwidths, in radians/sec. 

RBW Array of reciprocal bandwidths, in sec. 

-4 

OMM4 Array of o> 

OM3DD The factor C w AtuA0 in the numerator of e.g. (5) of the 

(X 

growth documentation. 

J X 16 

OM116 Array of w ’ , used in approximate computation of 

resonant wave number. See eg. (9) of §3.3.7. 

— 2 5 

OM25 Array of w , used in computing the exponent s in 

Miysuyasu's spreading function. See eg. (17) of §3.3.7, 

-4 -4 

DOMM4 Array of Used in cutting back high frequencies. 

Step 13 of §3.3.7. 

2 

OMK Array of wave numbers o> /g. 

3. 3. 4.6 XXGRID 

LAT coded latitude of current grid point. 

LONG coded longitude of current grid point. 

KPOINT 76+LAT+LONG. Used in addressing land-sea table. 

3 . 3 . 4 . 7 XXLAND 
LANSEA See §3.2. 

3 . 3 . 4 . 8 XXWIND 

WIND See DDl 0 in §3.3.1. 



26 


3. 3. 4.9 XXXH13 

JFRONT, HI 3 See DD14 in 53.3.3. 

3 . 3 . 5 . 0 XXPROP 

JPROP See account of JTAB4 in 53.1.5.7. 

PROP See account of TAB 4 in 53.1.5.7. 

3 . 3 . 5 . 1 XXPACK 

EEE, SPEC16 See DD13 in 53.3.3. 

3.3.6 Subprograms 

3.3.6. 1 MAIN. Declares common blocks and calls WWORK. 

3. 3.6.2 Subroutines called once only 

3. 3.6.2. 1 WWORK. Validates input, calls LODTAB, acquires an 

initial wind field, and sets up a do- loop to control 
growth, propagation, and output cycles. 

Called by: MAIN 

Calls: LODTAB, WWAX 1 , BUMP, QRWAX2 , and utilities. 

Method : 

1. Set default values of parameters (ISN 7 to 30). 

2. Read and echo namelist $WHAT (31-32). 

3. Validate input (33-165) . 

4. Acquire land-sea table (166-167). 

5. Acquire initial wind field (168-173). This wind 
field is used once only, in the first call to 
WWAX1 . See 53.3.7.1. 

6. Set ZMDI-I to the value corresponding to KSTEP1-1 , 
calling BUMP if KSTEP1 > 1 (174-179). 

7. Call LODTAB (180). 

8. Acquire the tables JPROP and PROP (181-184). N.B. 

In this SOWM the entire propagate table is treated 
as resident; in models written for smaller machines, 
the table was read in, repeatedly, one parallel 

at a time. 

9. Within do-loop in ISTEP, call WWAX1,BUMP, and 
QRWAX2 (185-190). 
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3. 3. 6. 2. 2 LODTAB. Computes the functions of direction and 
frequency contained in common block XPARAM. 

Called by: WWORK 

3. 3. 6. 3 Subroutines called at every time step . 

3. 3. 6. 3.1 BUMP. Updates the year, month, day and hour shown 
in ZMDH by adding the hours in DELTAH, and, if nec- 
essary, adjusts the result by adding 1 day and sub- 
tracting 24 hours (ISN7) ; adding 1 month and subtracting 
28, 29, 30 or 31 days (9-18); adding 1 year and 
subtracting 12 months (18, when MONTH = 12); sub- 
tracting 100 years (19-20) . 

Called by: WWORK 

3. 3. 6. 3. 2 WWAX1 . Reads wave spectra from DD17, applies half a 

time step (5400 seconds) of growth, and writes spectra 
to DDl 8 . The code at ISN28-35 copies spectra onto 
the east and west margins of the cut cylinder, where 
they are expected by the propagate routine. For the 
sequence of wind field acquisition in WWAX1 see §3.3.7. 1. 
Called by : WWORK 

Calls: UNPAQQ, CMPE27,PAQQ, and utilities. 


3. 3. 6. 3. 3 QRWAX2 . Reads wave spectra to DDl 8, unpacks them in- 
to array SPOLD, propagates to find the new spectra 
for one parallel and puts them in array SPNEW, applies 
half a time step of growth, writes the spectra at 
selected time steps to DDl 3, and writes the spectra 
at all time steps to DDl 7. The mathematical method 
is discussed in §3.3.8. The array MOOD, where 
MOOD(I) = MOD (1+1 , 3) +1 , is introduced to save computing. 
Called by : WWORK 

Calls: GRAB ,CMPE27 ,FAQQ and utilities. 

Method : 

1. Set up array (ISN7-10). 

2. At selected time steps, write ID to DDl 3 (11). 

3. Rewind DD17 and DD18 (13-14). 
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4. Acquire wind field (15). See 53.3.7.1. 

5. Zero the significant height field (16-18). 

6. Acquire two parallels of old spectra (19-20). 

7. In a do- loop on LAT: (21). 

8. If LAT < 72, acquire one parallel of old spectra 
(23-24) . 

9. Zero array of new spectra (25-27). 

10. In a do- loop on LONG: (28) . 

11. If the point is land, get out (29). 

12. Propagate. In order to follow the action at 
ISN32-42, it is advisable to compare each of the 
2+KPART+1 references to array JPROP and the KPART 
references to array PROP with the corresponding 
values printed in subroutine TRANSP (53.1.5.7.) 

13. Call CMPE27 to effect the second half step of 
growth (45) . 

14. Compute significant height (46-50) . Here endeth 
the loop on LONG. 

15. At selected time steps, write spectra to DD13 (55-56). 

16. At all time steps, write spectra to DD17 (57-58). 

Here endeth the loop on LAT. 

17. Write significant height field to DD14 (60). 

18. At selected time steps, write file mark on DD13 (61). 

19. Print significant height field (63) . 

3. 3. 6. 3. 4 CMPE27. Takes a spectrum and applies growth and 

dissipation for the time step DELT contained in com- 
mon block XXTIME. Note, that as used in this SOWM, 

DELT is half a propagation step, i.e. 5400 seconds. 

Called by: WWAX1 , QRWAX2 

Data: UXMULT(I) = u* , where U = (1+2.11) m/sec. 

U 1 9 

1 9 

UKAPPA (I) = U K , where U = (1+2.11) m/sec, 
19 1 19 

and k is defined in eg. (9) of 53.3.7. 

1 
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GAMMAK(I) » .4lK°' 28 f where k « .01*(I+1), 

, „ - 0.0528 

DELTAK(I) - .61543298k , where 

k * .01*(l+1) . 

These arrays are used to approximate the non-linear 
functions u*(U ig ), k (u*) , y(k), 6(k) by table look- 
up and interpolation: see XSN24-28, 159-160. 

The steps below are keyed to 13.3.7.2: 

1 . not specifically modeled 

2. (ISN16-17) 

3. (24,27) 

4. (25,28) 

5. not used 

6. (35-61) 

7. (62-63) 

8. (64-67) 

9. (72-112) 

10. (113) 

11 (114-115) 

12. (117-120) 

13. (122-128) 

14. (140-138) 

15. (139-145) 

16. (148-185) 

17. (187-211) 

18. (121, 212-216) 

3. 3.6.4 Packing and unpacking routines 

These subroutines were not part of the SOWM as first written: 
they were added so that DD17 and DD18 could each fit in 241 tracks. 
At such time as 478 tracks each are available for DD17 and DD18, 
some CPU time can be saved by reverting to subroutines WORK, 

WAX1, PRWAX2 . 


3. 3. 6. 4.1 PAQQ. Converts the spectrum at 76 grid points (85120 
bytes) into 42864 bytes and stores them in arrays EEE 
and SPEC16 in the format sketched for DD13 in §3.3.3. 
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Called by: WWAX1 , QRWAX2 

3. 3. 6. 4. 2 UNPAQQ. Reverses the operation performed by PAQQ, 

Called by: WWAX1 , GRAB 

3. 3.6. 4. 3 GRAB. Reads 85120 bytes from DD18, calls UNPAQQ to 
unpack the spectra into SPNEW, and then moves them 

to the addresses in SPOLD controlled by the parameter I. 

Called by: QRWAX2 

Calls: UNPAQQ and utilities 

3 . 3 . 6 . 5 Utility routines 

The purpose of these 6 subroutines is to segregate I/O and 
related statements, thus minimizing the work needed to replace, 
if desired, FORTRAN reads, writes, rewinds, and file marks by 
calls to non-standard routines. 

3. 3. 6. 5.1 BREAD2. Reads NPOINT 2-byte words from LU to ARRAY. 

3. 3. 6. 5. 2 BREAD4. Reads NPOINT 4-byte words from LU to ARRAY. 

3. 3. 6. 5. 3 BRITE2. Writes NPOINT 2-byte words to LU from ARRAY. 

3. 3. 6. 5. 4 BRITE4. Writes NPOINT 4-byte words to LU from ARRAY. 

3. 3. 6. 5. 5 BEND. Writes a file mark on LU. 

3. 3. 6. 5. 6 BREW. Rewinds LU. 
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3.3.7 The spectral growth and dissipation algorithm . 

The growth model used is the linear model 

dT/dt * A+BT (1) 

where T * SAcoA0 is the spectrum integrated over one of the 

(NDIRECxNFREQ) * 280 computational bands. Growth is controlled 

by not permitting T to grow above 

_ 2-4 4 4 4 4 

T = M(co,0){%ocg w [exp (-to /co ) -exp (-to /to ) ] '} 

0 0 02 01 /os 


where M(co,0) is an angular spreading function derived from 

Mitsuyasu et al, (1975)? the expression in {} is a Pierson- 

Moskowitz spectrum? < co^ are the limits of a frequency band, 

and u *= /(co to ) , As formulated, A/(AcoA0) is a function of w, 0, 
1 2 t 
and u* (the friction velocity) ; B is a function of the same 

arguments expressed as 


B = ioB*(om !t . cos 0 / g) ? 


(3) 


to is a function of U, the wind speed at 19.5 meters, 
o 

co » S^g/U; (4) 

0 

a is a function of e , the ratio of total energy to Pierson- 

*’** 1 *”* 2 4 

Moskowitz energy, e = E/(^a 0 g U ) ? M is a function of 0 6 

to/co . We consider each of A, B, co , a, M in turn. 

0 0 

A- term . Following Snyder 6 Cox (1966) and Inoue (1967, p. 

12) we write 

4 3 

C Ua-to AcoA0 

a * 

A = 

2 2 2 

[y+(k cos 0 - k) /y) [6+k sin 0/6) (5) 

where k = co /g, Y = MAX(.41 k ,.0027422197), 6 “ y /• 84 , 6 
k = co/u . The determination of u„ is the principal difficulty 
in applying (5) to wave growth. Inoue (1967, p. 22; u is a 
computation slip for u K75 ) took u c as the anemometer wind 
at 19.5 m. 
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According to Priestley (1965, p. 45) u_ is the wind at a 

C 

height z of 5,3 wavelengths, i.e. 

2uu_ 


z * 


2it 


5 . 3 k 5 . 3to 

To routinely solve the non-linear equation 


( 6 ) 


- * 2,5 u* log — ^ — 
* 5 . 3kz 

0 


for k as a function of w and u + is awkward, and the following 
device was adopted, based on locally replacing a logarithmic 
profile by a power- law profile; 

Take a reference height 

2 , 

z « 637. 39u li ,/g 


and compute 
> •* 

: 


v. tu 


2 -n/5 . 3z , 

Z - 1 ~ 1 

,00001525m sec u + +. 
2.5u*log(S/ Z ), 


-1 2 2 

001468m sec u^— . 0000371m, 

(7) 


(CU. 


For u + = 1 m/sec, corresponding to U g = 23.773 m/sec, this 

calculation yields 

z = 65.0 m 

k - .018238564 m” 

z = .00144615m 
0 

u - 26.783094 m/sec 

- 1 

a> * .48848518 sec or .0774483 hz. 

So far there is no approximation; to convert to other values 

of ui we introduce the scaling law 

- 7.25 , „ 

z/z - (u/u) . (8) 

Substituting (8) in (6) yields 

(u/u) ’ = (k/k) = (u/u) ((i)/u) ) 


„ -0.16 
u/u = ( oj/tt) ) 

- 1.16 

K = K (tt/co) = K (u*) XK (co) , 

1 2 


so that 
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where 

1,16 

K * KW 6 K •* U) 

1 2 


( 9 ) 


The constant C & was chosen by fitting the absolute power 
density scaling factor of Priestley (1965, p. 51) to u*, as- 
suming a u* law and roughness parameter z of .03 m 
corresponding to mowed grass) for the conditions of Priest- 
ley's experiment. The product of the fitted constant and all 

physical constants appearing in the relation between A and the 

- 9 2 

pressure spectrum is C a * 1.36x10 for A in m and t in sec. 
B-term . Introduce the non-dimensional phase speed 
t|j * wu* cos e/g } 

then B = wB* , where B* is a universal function of iji (Miles 
1959). Inoue (1967, p. 32) adopts the form 

-l z 2 ” 2 

B + = ( 2 it ) [.00139 exp{-7000 (i|i-. 031) } + .725i|i exp (-. 0004 <j, )). 

° ( 10 ) 

Two forms of B* have been implemented in subroutine CMPE26 : 

2 

B* = max (0 , p[i|i -.0004] ) , (11) 


where the choice B = .1155 yields the same asymptotic behavior 
as (9) for large i|»; 6 

B* = 0[ (^ 2 -.0004) + (503.3+{2042000+122040000^}{t|j-.031} 2 )" 1 ] (12) 

constructed as a rational approximation to (10). The form (12), 
with B = .1066, has been successfully used in duration-limited 
growth tests. 

Fierson-Moskowitz spectrum . The total energy in a "fully 
developed" spectrum is 


2 - 1-24 

%ag = ' l 4 aB g ,U , 


(13) 


where we adopt Pierson 6 Moskowitz' values a 
For immature spectra, Resio 6 Vincent (1977, 
the law 


2 - 4-23 

= .037 (Eg u* ) ’ . 



.0081, 8 a .74. 
19) propose 


et 


(14) 
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The exponent -.23 is a least-squares fit; Hasselmann et al. 

(1973, p. 37) found -.22; the parametric model of Hasselmann 

et al. (1976, p, 213) prescribes the exponent -.20. In order 

-4 

to preserve consistency with (13) we scale total energy by U : 

(15) 


-.2 3 -1 2-4 -.23 

a « a e ** a (4a Bg U E) 

0 0 0 


*"* 1 2 *“* 2 3 9 2 
with U in m sec and E in m , a « .0007293789E ’ U* 

Angular spreading . Mitsuyasu et al. (1975) obtain a cred- 
ible fit to the spreading function proposed by Longuet-Higgins 
et al. (1963) : 


M{«,0) = G(s) (k+h cos G) 


(16) 


where s is a function of 5 G(s) is a normalizing factor 

such that 

/M(cu,e)de » i. 

Mitsuyasu finds that s has a maximum at a)p PaJ< , ** . 8* 4 w , 6 that 

■ 2 . 5 

(17) 


s psak “ 11 - 5 ( Vak u 10 /s) 


where u is the wind speed at 10 meters. Thus, for the 
10 


- 2.5 


reference spectrum, 

Speak ’ 11 - 5 [( Vak / “ 0 > ( “ 0 u * /g)(u I() /u *n 
so that a value of u^/u,,, is needed to compute s pea j < . The 

range of u underlying (17) is 7 to 10 m/sec; we adopt 
l o 


u„ 


0.34 m/sec 


as a representative value of u # in that range. Then 


u 


l o 


2.5u*log(10/z ) , 


where z is computed from (7); and u /u* = 27.3484642 (so 
0 1 0 * 

that u = 9.2985 m/sec); whence s__ . = 15.00496. 

1 o peax 

For w 7 * , Mitsuyasu finds 

^ 5 —25 

s / s p e ak = mint( “ / Vak> ' '“/“peak' ' 
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but our Interpolation scheme for M gives difficulties for 
s < 1 and so we adopt the forms 

s ■ max 1 1 , 15.00496xmln{ (w/w penk ) * (w/w peak ) ’ )J 


(IB) 


and 


M(w,0) 


(1+cosG) s - 1 

. _ * 



[ (1+COS0) 


8 -1Jde 


(19) 


It was found that separating the algorithm into growth, cut- 
back and angular redistribution produced oscillations in the 
angular spectrum at high frequencies: accordingly, these oper- 
ations are interleaved inside an outer loop that cycles through 
frequency bands. 
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3. 3. 7.1 



The table below, covering the 12 growth cycles of 6 time 


steps, sufficiently illustrates the interaction of a 1^-hour 


growth step with wind 

fields 

archived 

at 2-hour 

intervals. 

ISTEP 

WAX 

t 

t+At 

t+%At 

wind field 

read Kind 

1 

1 

0.0 

1.5 

0,75 

1 

0.0 

1 

2 

1.5 

3.0 

2.25 

2 

2.0 

2 

1 

3.0 

4.5 

3.75 

3 

4.0 

2 

2 

4.5 

6.0 

5.25 

4 

6.0 

3 

1 

6.0 

7.5 

6.75 

- 

6.0 

3 

2 

7.5 

9.0 

8.25 

5 

8.0 

4 

1 

9.0 

10.5 

9.75 

6 

10.0 

4 

2 

10.5 

12,0 

11.25 

7 

12.0 

5 

1 

12.0 

13.5 

12.75 

- 

12.0 

5 

2 

13.5 

15.0 

14.25 

8 

14.0 

6 

1 

15.0 

16.5 

15.75 

9 

16.0 

.**■ 

o 

2 

16.5 

18.0 

17.25 

10 

18.0 


With the numbering here given, when prolonging a simulation, 

LSKIP * INT(3*KSTEP1/2)-1 


3. 3. 7. 2 


Sequence of Operations in growth anc 


jation algorithm: 


1. Read U and 0 W at grid point. 

2. If U < 3.11 m/sec (« 6.045 knots) exit. 

3. Compute u* from U by solving the equations 

u * » . 4U/log e (19.5/zo) / 

2 

z = o/u*+c u* +c , 

01 2 3 2-l -12 

where c = .00001525m sec , c = .00001468m sec , 
l 2 

c = -.0000371m, and c is chosen to agree with Garratt 
3 2 

(1977, p. 922) . 

4. Compute from u* by (9). 


5. Adjust e w for angle between true north and grid north. 
In the present SOWM these two directions are identical, 
and no adjustment is made. 



37 


6 . 
7, 


Compute cos(cT-e w ) for midpoint ? of each angular hand used. 
Compute * . 74^g/U and *» . 


8. Compute the total spectral energy E (all frequencies, 
all directions) . 

9. Compute the virtual reference spectrum a by (2) , 


(17), (18), for all frequencies and all downwind directions. 
[N.B, No specific tables of upwind and downwind directions 
are kept? the sense of a direction is determined by 
examining the sign of cos (?-e w ) . 1 

10. The remaining steps (11-19) are enclosed in a loop that 
cycles through frequencies, beginning at the lowest. 

1 1 . Compute a by (15). 

12. Compute the one dimensional spectral component 

T (w) * ]> T (w , 6) (all directions), 

- 2-4 - 4 \ 

13. if T(») < ^ag (to < -to n '* cycle through directions to 

k fit 

; T new (w,0) - T(w,fl) xka^ (w^-oT^/Tdu) . 

14. Dissipate opposing bands . Cycle through pairs of direc- 
tions {0,0+rr}: A « | T (e, 0) -T (w , 0 +it) | /{T (w , 0) +T [in , 0+ir ) ) ; 
T new ( w,0) » AxT(w, 0), T new (w,0+ir) » A*T (w , 6+ir) . 

15. Compute the downwind energy 


T . « £ T(w,e) + - £ , T(&r,0). 

|6-0 w | < }£tt 2 |O,0 w l= J 5ir 


2-4 4-4 4-4 

If T. < hag ot texp(-w w )~exp(~cu <u )] (underdeveloped 
« 0 0 1 0 2 

band) , go to step 16; otherwise to step 17. [In this 
model, at any time step, a frequency band undergoes 
growth or angular redistribution, but not both.] 

16. Grow . Cycle through downwind directions: if 

T(«,6) > T (w,6), or B*(w,0) < 0, 

do nothing, otherwise 
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T new (w,e) - MINIT (w,e) ,{T(w,0)e BAt +Alf 1 (e BAt -1)}] . 
Go to step 18. 

17. Redistribute energy over angles. Por each frequency 
band, cycle through downwind directions and compute 

x(w,e) « MiN[T(u,e)-T (uT,e),oj, 

0 

Y (<o,0) «» MINI? (w,0)~T(w,0) ,0] . 

0 

Sum X and Y to obtain X d and Y^. Compute 

5 = X d /T d and n « V Y d T d* 

Then 

T new (u,0) = T(w,0) - 5X(w,0) + nY(w,0). 

18. Adjust the total energy E: 


E « E , , - T , j ( w ) + T ( a) ) . 

new old old v ' new v ' 


3,3.8 The spectral propagation algorithm 

The physics of wave propagation in deep water is well under- 
stood as a result of the work of Barber 6 Ursell (1948), Groves 
6 Melcer (1961) and Snodgrass et al. (1966). Each component in 
the two-dimensional spectrum travels along a great circle in its 
direction at the deep-water group velocity appropriate to its 
frequency. The development of an efficient and accurate computer 
algorithm to simulate this process has been particularly difficult. 

Baer (1962) studied the simple first-order finite-difference 
analogue to the convective term (velocity-gradient technique) 
and concluded that it was inadequate if the quasi-discontinuous 
spatial distribution of wave energy was to be preserved. Such 
a scheme has nevertheless been used, for example by Barnett 
(1968). More recently Ewing (1971) has applied a fourth-order 
convective difference scheme, which he showed to be considerably 
more conservative than first-order schemes. A recent paper by 
Brian Golding (1979) indicates that the latest practice at the 
U.K. Meteorological Office, Bracknell, is to use a second-order 
Lax-Wendroff scheme with a polar stereographic grid. 
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The so-called jump technique was proposed by Baer (1962) 
to overcome the smoothing effects of the velocity-gradient 
technique. In the jump technique, wave energy is considered to 
be discon tinuously distributed on the grid system; and spectral 
components are simply translated (jumped) to adjacent grid 
points after a sufficient number of time steps has elapsed to 
account for the quotient of the spacing of grid points by the 
propagation velocity. The FNWC operational Mediterranean 
spectral wave model described by Lazanoff, Stevenson $ Cardone 
(1973) employs the jump technique. Uji 6 Isozaki (1972) have 
developed a more complicated version of the jump technique 
whereby lateral spreading and longitudinal dispersion associated 
with discrete directional spectral components are simulated. 

The postulate central to the development of all schemes 
described above is the faithful propagation of monochromatic 
waves. However, waves on the ocean are not monochromatic; even 

a narrow-band swell has finite bandwidth. The error incident 

k 

to a crude gradient scheme xs of order N , where N is the number 
of time steps; the error inherent in neglecting bandwidth is of 
order N, and is the dominant error for long propagation distances. 

A conservation property considered essential has not been 
rigorously verified for schemes described above, and is suscep- 
tible of a priori calculation only on a flat earth: indeed, the 
numerical experiments of Uji 6 Isozaki (1972) did not extend 
beyond consideration of a flat earth. The property, which may 
be called integral conservatio n, states that, for any distribution 
of energy removed from land and from the edges of the map, the 
total energy integrated over all grid points shall be identically 
equal before and after a time step. The scheme submitted 
achieves integral conservation by the exclusive use of downstream 
interpolation. 

After the propagation matrix has been calculated as sketched in 
S 3 . 1 . 6 , and read by the SOWM in the format stated in § 3 . 1 . 5 . 7 , the 
propagation proper reduces to KPART address calculations to check 
whether the transmitting point is land or sea, another KPART 
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address calculations to locate the transmitting spectral 
bins in the array SPOLD, and KPART multiplications and additions 
to compute the variance in the receiving spectral bin in SPNEW. 

See 03.3.6.3.3, step 12. 

All propagation distances are less than 120 miles, guar- 
anteeing stability for north-south propagation, The occasional 
jumps into triangles XIX and IV at high latitudes and low fre- 
quencies have not been extensively investigated for stability 
with this SOWM? but similar constructions in SOWMs for limited 
basins have caused no trouble. In a globe completely covered 
with water, the approximation adopted for convergence of merid'- 
ians would lead to a slowly amplifying band of wash at the 
equator in directions 81°, 99°, 261°, 279°; in the model as 
submitted, the growth of this wash is blocked by sinks in 
Africa and South America. 

3 • 4 Readback programs 

3.4.1 DUMP13. Reads spectra from tape 13 in the format in 

which they were written by PRWAX2, and prints the 
spectra, with marginal totals, for a selected list 
(not more than 100) of grid points. 

Input parameters (in calling sequence) : KPOINT is the 
number of points to be dumped. KSTEP1 , KSTEP2, KSTEP3 
have the same meaning as in the subroutine WORK 
§ 3 , 3 . 6 . 2 . 1 ) ; they must have the same values they had 
in namelist $WHAT. 

Input parameters (in namelist $ WHERE) : LPOINT is the 
list of points to be dumped. The points must be in 
ascending order; no check is made for land or for 
LLONG = 1, 2, 75, 76, but any point with those properties 
will produce a page of meaningless output. 

3.4.1 DUMQ13. Reads spectra from tape 13 in the format in 

which they were written by QRWAX2 (§3. 3. 6. 3. 3), and 
prints the spectra, with marginal totals, for a se- 
lected list of grid points. Input parameters and usage 
are the same as for DUMP 1 3 . DUMQ13 does not use sub- 
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routine UNPAQQ (§3. 3. 6. 4. 2); only spectra at points 
to be printed are converted from integer to real. The 
mathematical method is the same as in UNPAQQ , viz 
reversing the steps, taken in PAQQ . 

3.4.3 DUMP 1 4 , Reads significant height fields from DD14 

in the format in which they were written by PRWAX2 
or QRWAX2 , and prints a table of coded latitude, 
coded longitude, and significant height. 



42 


4. APPLICATION OP MODELING STUDIES TO SEASAT DATA 
4.1 Assimilation of Scatterometer Derived Surface Winds 


It is a characteristic of all operational numerical weath- 
er prediction models and of most general circulation models, 
including the GMSF GCM, that planetary boundary layer (PEL) 
exchange processes are treated quite parametrically. That is 
fluxes of momentum, heat and moisture across the PBL are param- 
eterized in terms of the wind, temperature and humidity at the 
lowest active level of the model which is usually near the top 
of the PBL, and surface temperature, roughness and moisture 
characteristics (boundary conditions) . In most models, there- 
fore, the near surface wind, which is the quantity sensed over 
the ocean by the SEASAT SASS, is not specified or forecast 
explicitly. (In the GMSF GCM, a surface wind is computed in- 
ternally through a simple downward extrapolation of the wind 
from higher levels.) Initialization in such models is concerned 
primarily with the mass field and the windfield above the PBL. 

The simulation studies described above avoided the problems 
of assimilation of surface wind data into models of the above 
type by synthesizing the simulated scatterometer winds at the 
lowest active level of the GMSF GCM, level 9, which is situated 
near the 950 mb level. The simulation study therefore implies 
perfect boundary layer physics, since no errors have been intro- 
duced in the simulation to account for the extrapolation of 
winds from the surface to level 9. Another simplification made 
in the simulation was the avoidance of the "aliasing problem", 
which derives from the fact that the SASS wind algorithm operates 
on colocated pairs of measured radar backscatter cross-section 
measurements to return up to four possible wind vector solutions 
at each cell on the sea surface sensed. While the wind speed 
differences between the wind vector solutions are typically very 
small, the wind directions differ significantly. In the four 
alias case, the directions may be nearly 90° apart. In the three 
alias case, two of the three vectors may be only 30° apart, while 
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the two alias solutions are usually 180° degrees apart. 

Experiments with real SEASAT-A wind data will therefore 
require some modifications of the procedures used in the simula- 
tion study, As a part of this study, a preliminary and rather 
simple^modification to the GMSP GCM was designed, coded and 
tested , The modification was intended to satisfy two objectives: 
(1) provide a reasonable specification of the marine surface 
wind direction at each GCM grid point, each time step, to be used 
as a first guess in a SASS alias removal scheme (the modification 
as presently coded simply picks the SASS alias closest to the 
first guess) ; (2) allow the assimilated SASS surface wind speed 
and direction data to be extrapolated to the lowest prognostic 
level of the GCM. 

The procedure devised to accomplish the objectives utilized 
a simplified version of the PBL model developed by Hoffert and 
Sud (1976). That model is used to diagnose the surface wind, 
defined as the effective neutral stability wind at the 19.5 meter 
height, from level 9 and surface information. The method accounts 
for baroclinicity in an approximate way by superimposing baroclinic 
turning of the wind in the PBL upon frictional turning in com- 
puting the surface wind direction. 

The basic procedure for specification of the surface wind 
may be outlined as follows: 

1 . At each ocean GCM grid point assemble the 
following quantities : level 9 wind components 
ug , v 9 ; level 8 wind components Ug , vg ; 

level 9 potential temperature 0 g , and humidity 
qg; sea surface temperature, 0 g . 

2. Compute a surface wind direction from level 9 
and level 8 wind components using the vertical 
extrapolation scheme presently used in the GCM. 

The difference between the extrapolated surface 
wind direction and the level 9 wind is considered 


Oceanweather Inc. contributed primarily to the development of a 
functional description, and to evaluation of simple tests; coding 
was done by M. Helfand of GMSF. 
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3. 


the approximate barociinic component, Ac* b , of 
the turning of the wind in the PBL. 

Compute the magnitude of the friction velocity, 
u + from the simplified Hoffert-Sud parameterization; 
that is, 

i—— 2 2 L 

U * \ C d X (U9 + v 9 ^ 

where the PBL bulk drag coefficient C d is com- 
puted as a function of the dimensional ratios 
h/zo end R b , where h is boundary layer depth, 
is the roughness parameter and R b is 

R b = gh (0 g-0 s ) 

2 2 ” 

0 s (u 9 +v 9 ) 


The calculation of proceeds in two stages. 
First, the neutral drag coefficient C^ n is 
obtained by iteration between 

c dn = - 15 /(log e h / z o) 2 


and 


z 

z 0 « A/u* + Bu* + C 


where A, B, C are taken from Cardone (1969). 

The first guess of z 9 is taken as .0025 cm and 
h is taken to be the height of level 9. The 
drag coefficient C^ n then is modified to account 
for stability through functions, Fj , which are 
curve fits derived by Hoffert-Sud from their 
numerical PBL solutions : 

C d s C dn * Fl (R b' h/zo) 

A separate set of curve fits, F^ , provide the 
frictional turning of the wind Aa f between 
level 9 and the surface: 

A a f = F 2 (H^, h/zo) 

Given Aa fa and Act f , the surface wind and stress 
direction are obtained directly. The effective 
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19.5 meter wind speed is defined simply by surface 
layer theory 

I V 1 9 . 5 ! ** 2.5u + log e (1950/z 0 ) 

The assimilation of SASS winds follows a slightly different 
procedure, as winds at the surface must be extrapolated upward 
to level 9. The assimilation proceeds in the following steps: 

1 . From the surface windfield computed as described 
above, interpolate surface wind directions to the 
locations of the SASS measurement cells, and 
where aliased solutions exist, choose the solu- 
tion closest to the interpolated direction. 

2. Since the SASS solution may be described 
equivalently in terms of' I V 9 . 5 I = 2 . Su^log (1 950/z q) , 
or u*, compute 

R b * * gh( 0 9 -e s ) 

and 2 

Z 0 = A/u* 4- BU* + C 

3. Employ the simplified Hoffert-Sud model expressed 
in terms of F 3 and F 4 to compute the magnitude 

of the level 9 wind from 

|V 9 | = u* 2 /C d 

where 

C d = F 3 (R b *, h/z 0 ) 

4. Compute the frictional turning &a^ from 

Aa f = F 4 (R fa * , h/zo) 

5. Compute the level 9 wind direction from the surface 
direction, the computed frictional turning Aa f , 

and the baroclinic turning Aa b of the guess field. 

The modifications described above have been tested mainly 
to assume error free code. No SEASAT simulations were performed 
with the assimilation algorithm. 
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4 . 2 Impact of SEASAT-A data on large scale weather forecasting 

The SEASAT-A data processing activity has progressed to the 
point where global SASS data sets will soon be produced, archived 
and disseminated. The study of weather patterns over the Northern 
Hemisphere by the SEASAT Simulations Studies Steering Group, 

(see Appendix A) has shown that the last month of the mission 
provides suitable conditions for the assessment of the impact on 
forecasts of data from SEASAT-A. An archive of all data necessary 
to perform the impact tests is being assembled at GMSF. The 
conventional analyses and data sources have already been assembled. 
The SEASAT-A data archive will consist initially of SASS data 
in terms of sensor data (backscatter coefficients) and geophysical 
data (wind vector aliases) . 

The first test of data impact, therefore, will be a straight- 
forward application of the methodology of the simulation studies. 
Objective analyses of the conventional data sets will yield an 
estimate of Nature. Control analyses and forecasts will be pro- 
duced exactly as in the simulation experiments . A SEASAT 
assimilation run will repeat the control assimilation and addi- 
tionally include the asynoptic assimilation of data from SEASAT-A. 

As described above, the GCM model winds, interpolated to the mean 
cell location, will be extrapolated to the surface and used to 
remove the aliases. The satellite wind vectors will be extrapolated 
back to the lowest active level of the GCM and assimilated by the 
SCM as in the simulations. Forecasts from the SEASAT assimilation 
run will be compared with control forecasts, with objective 
analyses produced by the GMSF analysis package, and with an ob- 
jective analysis produced by an analysis system independent of 
the GMSF GCM. If SMMR measurements are available they may be 
used to correct SASS data for attenuation by clouds and rain 
where possible; the SASS data will be rejected where such correc- 
tion is requisite but not possible. 

It may happen, contrary to expectation, that the first tests 
with real data from SEASAT will confirm the results of the 
simulation studies and show large beneficial impact on short- 
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range numerical weather forecasts. Rather, we expect the first 
tests to uncover problem areas requiring the development of more 
intricate experimental techniques. We anticipate at least one 
of the following scenarios, not to mention some combination of 
them: [1] the scatterometer data are not as accurate a measure 
of surface wind as expected; [2] the scatterometer is apparently 
a good measure of surface stress, but the parametrization of the 
planetary boundary layer in the GCM is too crude to utilize the 
data; [3] the use of the GCM forecast in the retrieval process 
results in the rejection of valid wind data from SEASAT-A; [4] 
the wind data are duly assimilated, but their effect on forecasts 
is not significant. 

If an initial test should show large beneficial impact, the 
program of experiments can proceed to develop optimal assimilation 
methods incorporating, in addition to SEASAT, other remotely 
sensed data such as indirect soundings and winds from cloud 
trackings; and thence to design an operational satellite system. 
Problem [1] would require more sophisticated preprocessing 
techniques, based on geophysical reduction algorithms more com- 
plicated than those now before the SASS team. For example, back- 
scatter may depend jointly on wind stress and large-scale surface 
roughness; to retrieve surface stress would require that a coupled 
GCM-SOWM be used in the assimilation. Problem [2] would place 
the emphasis on model development; as part of the simulation 
studies program, a hierarchy of modifications to the parameter- 
ization of the planetary boundary layer at GMSF has been formu- 
lated; while only the simplest modifications are being implemented, 
more complicated models can be made available as needed. Problem 
[3] would require more complicated assimilation procedures. One 
possibility is to increase the assimilation period to 3 or even 
6 hours; so that several orbits of SEASAT-A can be processed 
jointly with ancillary satellite data (e.g. wind vectors from 
low-level cloud tracking) , conventional ship and buoy reports, 
and model forecast information to provide an improved guess field 
for alias removal. A man-machine mix analogous to the "special 
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effort" procedures being implemented for the global weather experi- 
ment is a conceivable operation at this step, but obviously it 
should be implemented only if found to be necessary. Problem (4) 
would represent a confirmation of the earliest suspicions of the 
Simulation Studies Steering Committee/ viz; that even accurate 
wind data should not significantly affect numerical forecasts be- 
cause of the slight contribution of the surface layer to the 
overall energetics of synoptic-scale systems and because/ in 
accordance with geostrophic adjustment theory, the model will 
reject the wind data dynamically. If so, surface winds from 
SEASAT would have to be combined with other data sets and, 
possibly external to the assimilation step, be extended to influence 
other variables, such as the pressure field. In any case, the 
joint assimilation of SEASAT winds and sounder information should 
be tried; a synergistic effect between the two kinds of data was 
apparent in an early simulation study. 

4 . 3 Utilization of Wave Data from SBASAT-A in GSOWM. 

There is already considerable evidence to support the validity 
of the wave heights to be retrieved from SEASAT-A. The concept 
of obtaining RMS wave height by analysing the shape of the return 
pulse to a radar altimeter has been proven? GEOS- 3 has been 
routinely returning wave heights so retrieved for several years. 

The improved radar altimeter on SEASAT-A is likely to attain its 
objective of accuracy of +0.5 m or +10$ in significant wave height. 
The wave height data are not now used operationally, though they 
are of potential benefit. Pierson & Salfi (1978) compared signi- 
ficant wave heights measured by GEOS- 3 with operational northern 
hemisphere wave analyses produced at the U.S. Navy Fleet Numerical 
Weather Center (FNWC) for 44 orbit segments obtained during 1975 
and 1976. The comparison revealed a significant bias in the wave 
analyses and occasional large differences over some orbit segments; 
the latter were attributed to poor specifications of windfield 
input to the wave model that generated the analyses. The wave 
heights from GEOS- 3 were not compared with forecasts. 

Hie basis of wave analysis and prediction at FNWC is a numerical 
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spectral wave specification model operated in a hindcast-f ©recast 
cycle twice daily, using windfield input data exclusively; in 
wave forecasting, unlike numerical weather prediction, initial 
wave conditions are not specified from wave data, which have 
always been too scarce to make such a specification profitable, 

GEOS- 3 and now SEAS AT- A, however, have altered the global wave 
measurement data base dramatically. As part of this study, a 
spectral wave specification model was adapted to the GMSF computer, 
to allow the development of methods for assimilating wave heights 
from SEASAT-A into a wave hindcast-f orecast cycle based on that 
model. 

Contemporary wave specification models have been reviewed 
recently by Cardone 5 Ross (1979). All such models are based upon 
the spectral energy balance equation, usually applied in its 
simplest form; that is, to surface gravity waves assumed to 
propagate through water of infinite depth that is otherwise at 
rest. In this form, the equation may be written 

I t E(f,0,x,t) + £ (f,e)xAE(f,e,x,t) 55 s(f,0,x,t) (20) 

where E is the energy density of the wave field described as a 
function of frequency f, direction of propagation 6, position x, 
and time t; is the deep-water group velocity vector; and S, 
the source function, represents all physical processes that 
transfer energy to or from the spectrum. Models of the so-called 
discrete type represent the spectrum E in terms of a number of 
spectral components (bands) of finite width, and successively 
simulate wave propagation (the homogenous part of (20)) and local 
energy transfers (the r.h.s. of (20)) in a series of discrete time 
steps on a grid overlaid on the basin of interest. Both rectangular 
and triangular grids have been used. In the model adapted to the 
GMSF, the spectrum is resolved into 192 components (16 directions 
x 12 frequencies), and the scheme reported by Greenwood 6 Cardone 
(1977) is used to propagate spectral components along archs of 
great circles and to refer the resulting propagation pattern to 
the GMSF GCM spherical grid. The source function S accounts 
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explicitly for energy transfers from the wind and implicitly for 
non-linear wave-wave interactions (Cardone et al. 1976) , A dis- 
crete model is run operationally at FNWC to provide 72-hour fore- 
casts for basins in the northern hemisphere/ twice daily. Initial 
conditions E(f,0,x,to) are calculated in a hlndcast procedure 
using the same model and surface windflelds produced by objec- 
tive analysis of ship reports of wind and sea-level pressure. 

An alternate model context for wave prediction, proposed 
recently by Hasselmann et al. (1976), is based on a parametric 
representation of the spectrum 

E (f , 0 ) * E(a^, A*1,...,n) 

involving a relatively few parameters a. . The JONSWAP spectrum, 
much used in this context, specifies the five parameters f ffl , a, 

Y, o a , and o b ? where f m is the frequency of the spectral peak, 
a is the "equilibrium range" constant, f is a peak enhancement 
factor, and a & and o b specify the width of the left and right sides 
of the spectral peak. The directional spread of energy is taken 
as symmetrical about the local wind direction, The parametric 
models are implemented by projecting the energy balance equation 
(20) including the source terms, onto a set of prognostic equations 
for the parameters of (21). Hasselmann et al. propose that for 
wind-generated seas only one or two parameters are necessary: their 
one-parameter model, for example, specifies a quasi-equilibrium 
relation between a and the non-dimensional peak frequency 
v ~ £ m u io/9/ where Ujq is the wind speed 10 m above the sea surface 
More complicated wave regimes require a hybrid model combining a 
parametric model of local sea with the propagation of swell in 
bands , 

While discrete models have been fairly extensively applied, 
calibrated, and verified, there have been few studies comparing 
discrete spectral models with hybrid-parametric models. It appears 
however, that some of the concepts advanced in the development 
of parametric models are suited to the development of techniques 
for assimilating remotely sensed wave data into a wave hindcast- 
forecast system. Significant wave height, the characteristic of 
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the sea surface measured by the radar altimeter on SEASAT-A, is 
an Integral property of the spectrum and is therefore more appro- 
priately assimilated in the parametric domain. It will therefore 
be necessary to incorporate some elements of parametric models 
into the GMSF SOWM for the purpose of assimilating significant 
wave height. The method will require, first, modification of the 
model so that at any desired grid point and time step the discrete 
spectral matrix can be partitioned into swell and wind sea? and 
the latter represented parametrically, in terms of total energy, 
peak frequency, and mean direction of wave travel. Then an algo- 
rithm may be developed that accepts updated values of total energy, 
wind speed, and wind direction, and adjusts the remaining parameters. 
In developing that algorithm by transforming the prognostic 
equations of Hasselmann et al. , it appears that a more convenient 
pivotal quantity car be chosen than the non-dimensional peak 
frequency v? Resio 6 Vincent (1977) show that the non-dimensional 
total energy e, which in simple fetch-limited growth models is 
linearly related to the non-dimensional fetch, is a useful 
substitute for the latter in more complicated regimes. 

The methods may be tested as follows. Given a run that 
assimilated real winds from SEASAT-A, use the surface winds speci- 
fied therefrom at 3-hour intervals to drive the GMSF global SOWM 
until the hindcast wave field is independent of the initial con- 
ditions in the SOWM (usually several days) . Then begin to assim- 
ilate significant wave heights from SEASAT-A. At grid points near 
the altimeter swath, partition the spectra into swell and para- 
metric sea, deflate the measured significant height by the swell, 
and update the SOWM sea on the basis of the sea inferred from the 
measurement. Then, by integrating the updated parametric sea 
over the 192 frequency-direction bands, reconstruct the directional 
spectrum at the grid point and add the swell back in. A method 
of successive corrections in parameter space may be applied to 
blend the corrections into grid points adjacent to the subsatellite 
arcs treated at the assimilation step. The method assumes that 
the difference between modelled and measured wave height at a grid 
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point is primarily associated with the wind sea and therefore 
with local and recent errors in the windfield. This assumption 
is based on the intuitive notion that the continuous insertion 
of global wind and wave height data, if successful, would cor- 
rect most errors in spectra in the zone of generation before they 
could be propagated any great distance. The assumption, however, 
needs to be tested. 
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5. SUMMARY AND CONCLUSIONS 

Within the next year, high resolution global data sets con- 
sisting of measurements of marine surface wind, stress and signi- 
ficant wave height derived from SEASAT-A microwave sensors will 
be processed and made available for evaluation. Those data sets 
will be truly unique in providing global coverage of air-sea 
interfacial parameters that are at best poorly sensed from sparsely 
distributed ground based systems. Of great significance is the 
fact that the scale of motion sampled in the satellite resolution 
cell matches closely the synoptic scale and therefore is of great 
potential value in the initialization of large scale numerical 
weather forecasting and ocean sea state prediction models. Con- 
ventional anemometer measurements for example, do not adequately 
separate the turbulent and synoptic scale because of the short 
averaging interval characteristic of even the best surface based 
systems. 

Preliminary evaluation of the performance of the SEASAT-A 
sensors of interest in this study (SASS ,ALT ,SMMR) suggests strongly 
that the design goals for measurement accuracy have been met. It 
remains to be demonstrated however, that skill in operational 
weather and sea state analyses and forecasts would increase if 
such data were available routinely from an operational SEASAT 
type satellite system. 

The studies reported here provide a firm basis for the conduct 
of an experimental program to assess the potential of real SEASAT 
data. Three major accomplishments are reported. First, as part 
of a collaborative program with the GMSF, a series of observing 
system simulation experiments has been performed to assess the 
potential impact of marine surface wind data on numerical weather 
prediction. The experiments simulated the time continuous 
assimilation of remotely sensed surface wind data following a 
SEASAT-A orbit. When error-free winds were assimilated using a 
localized successive correction method, substantial impacts in 
simulated 72 hour forecasts of surface pressure were found over 
both land and ocean extratropical regions. The effect of nominal 
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SASS errors on the impacts were found to be small. The simulation 
studies suggesi', that marine surface wind data if accurate enough, 
can have a beneficial effect on numerical forecasts comparable 
to or even greater than the effect of accurate remotely sensed 
sounder data. 

The second accomplishment was the development and implementation 
on the GMSF computer facility of a global spectral ocean surface 
wave specification model of contemporary formulation. The program 
models the generation, propagation and dissipation of surface 
gravity waves on a globe defined by the GMSF fine grid. Thus 
the model may be run in tandem with the GMSF GCM in either of 
three modes. In one mode, the model can be driven simply by 
analysis and forecast windfields specified by the GCM. In such 
a mode, the model serves to translate satellite induced impacts 
in surface windfield analyses and forecasts to surface wave 
analyses and forecasts. In another mode, the wave model itself 
might be used to assimilate remotely sensed significant wave 
height to reduce errors in initial wave spectra. Finally, should 
refined SASS geophysical retrieval algorithms require the knowledge 
of sea state for the optimum recovery of the surface wind and 
stress field, the GCM and GSOWM may be coupled to allow the 
retrieval and assimilation of surface winds of maximum attainable 
accuracy from scatterometer measurements. 

Finally, this study included the design of algorithms 
for the assimilation of SASS surface data into the GMSF GCM 
and for the utilization of real SASS wind data and ALT wave height 
data into a coupled GCM-SOWM. The former algorithm has been 
implemented in an experimental version of the GMSF GCM to be used 
for the first real data tests. The latter algorithm is largely 
in the conceptual stage. 
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Appendix A: Papers, presentations, meeting participation supported 

under this contract. 

Papers 

1 , Realistic simulations of the global observing system and 
of SEASAT-A marine wind data. M, Cane, V.J. Cardone, 

I. Halberstam, M. Halem, J, Ulrich. Submitted to 
Mon. Wea. Rev. 

2. Observing system simulation and potential impact of marine 
surface wind data on numerical weather prediction. Same 
authors as above. This paper is a substantial revision 

of the above paper with considerable material added. 
Submitted to Mon. Wea. Rev. in April 1979 as a revision 
of 1. 

Presentations 

1 . Realistic simulations of the global observation system 
and of SEASAT-A marine wind data. M. Cane and V. Cardone. 
Presented by M. Cane at the 1978 Pall Meeting of the 
American Geophysical Union, San Francisco, Calif. 

(Abstract attached) . 

2. Results of SEASAT simulation studies. V. Cardone. In- 
formal presentation to the interagency NOSS Advisory 
Panel, July, 1978, Jet Propulsion Laboratory, Pasadena, 
Calif . 

3. The contribution of SEASAT satellite capabilities to 
marine forecasting. V. Cardone. Presented to Inter- 
national Maritime Weather Conference and Exhibit, 

February 7, 1979, New York. 

4. Results of GLAS SEASAT Simulation Studies. V. Cardone. 

4 

Presented to S G SEASAT Planning Meeting, June 6,7, 1979 
Atmospheric Environmental Service, Downsview, Ontario 



Meetings 

1. SASS Team Meeting, Langley Research Center, 

August 15-16, 1978 

4 

2. S G SEASAT Planning Meeting, June 6-7, 1979, AES 
Downs view, Ontario 

3. Travel to GMSF, Greenbelt, Md. by Oceanweather staff. 

1978: April 13-14, May 24-25, June 26-27, August 2, 

August 28, September 26-27, November 14-15, December 20-21. 
1979: April 5-6, April 24-25, July 5-6, August 28. 


REALISTIC SIMULATIONS OF THE GLOBAL OBSERVING 
SYSTEM AND OF SEASAT-A MARINE WIND DATA 

Mark A. Cane (NASA/Goddard Space Flight Center, 

Greenbelt, Md. 20771) 

Vincent J. Cardone (Oceanweather Inc., White 

Plains, N.Y. 10601) 

A 30-day history run made with the Goddard 
Modelling and Simulation Facility General Circ- 
ulation Model (GCM) was used to fabricate simu- 
lated observations at the times and locations 
of the conventional surface, radiosonde and 
ship reports actually received during Feb. 1976. 
The fabricated observations, suitably degraded 
for instrument and sampling errors, were then 
used to create analysed fields on the GCM grid 
in an analysis-forecast cycle like those in use 
at major meteorological centers. The control 
fields so produced are much more representative 
of actual analyses than those produced by per- 
turbation of initial states with random errors. 
Significantly, the forecast error growth in 5 
simulated 72-hour forecasts from the control 
states is similar to that found in operational 
numerical forecasts. Further experiments simu- 
lated the addition of surface winds derived 
from SEASAT-A to the control run. These winds 
were fabricated directly for the GCM grid 
points intercepted by a SEASAT-A scatterometer 
swath, and were assumed to be error-free and 
representative of the lowest active level of 
the GCM. Two different asynoptici assimilation 
methods were used: direct insertion (DIM) and 
successive correction (SCM). Assimilation of 
winds by DIM was found to produce slight im- 
provements in simulated analyses and forecasts. 
The SCM experiments, however, resulted in 40% 
reductions of error in the specification and 
72-hour forecast of sea-level pressure and low- 
and mid- tropospheric winds. Preliminary results 
will be presented from experiments in progress, 
which include the simulation of errors expected 
in scatterometer-derived winds; and, as avail- 
able, from experiments using real SEASAT-A 
global wind data sets. 
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PROGRAM PRELIMS/MAIN 


COMMON 

S/ VP ARAM/ LLAT.LLONG.LPNT1 ,LPNT2,FREG(14) .FREQ1 t FREQ 2* SPACE ( 2 ) 
S * WEDGE * D 1REC( 20 ) »DELT 
S/VSHORT/ TLAT(3 .73) , TLCNG(2.76> 

S/YYOLTB/ JTABLE. TABLE! 4,20,73) 

S/YYTRIG/ DLAT.COSS.D ANGLE .CLAT . TRAVEL. RAD.CTR »STR, 

* 0MIN.O17,RAT(3.73) 

IMEGER42 JTABLE ( 16 . 20 .73 ) 

REAL* 8 D LAT ( 73 ) ,COSS (20 ). DANGLE ( 20 ) .CLAT ( 2 .73 ) ,TR AVEL ( 1 A ) , 

* CTR ( 1 A ) • STR (14) »RAD»0MIN.Q17 
DATA L13.L14/13.14/ 

FREQ 1 = .04 
FREQ2 = .20 
*EDGE = S. 

DELT » 3. 

CALL GGRID (L)3) 

C ENDFILE LI 3 

CALL BANDS 
CALL TRIG 
DO 10 IA - 1.14 
CALL OUTBND (I A) 

CALL INB.ND (L14+IA) 

10 CONTINUE 

CALL TRANSP (L14.L13) 

ENDFILE L 1 3 

STOP 

END 


U1 uu 


D2 


SU3RCLTINE GGRID <LU» 

COMMON 

S/YPARAM/ LLAT.LLONG.LPNTl ,LPNT2 .FRECC 14 > , FREGI ,FREQ2» SPACE<2 ) 
* , WEDGE ,D IRECI 20 ) ,DELT 

*/ YSHORT/ TLAT(3.73) ,7LQNG<2,76) 

S// JLA7 • JLCNG » ALAT( 76,73) *AL0NG(76 ,73 ) 

NAVEL. I ST 

S /LI ST 0 2/ LLAT.LLONG.LPNTl .LPNT2 .SPACE 
IN7EGERA2 JLAT<76,73 ) , JLONG ( 76 ,73 ) 

DATA LCUT/6/ 

REMIND LU 
LLAT = 73 
LLCNG = 76 

LFNTl = LLAT*(LL0NG--4) 

LPNT2 = LLATALLONG 
SPACE ( 1 ) = .999848 
SPACE (2) = . 325568 
DO 27 IA ” I .2 

SPACE ( I A ) = SQRTO 20.**2+( 1 50 . ASP ACE ( I A) )**2) 

27 CCNTINLE 

WRITE (LOUT ,LI ST 02) 

DO 22 IA = I ,73 

TLAT ( 1 , I A) = 2*IA-74 

TLAT (2 , 1 A) = 300.*C0S(.0174533*TLAT<1,IA>) 

TLAT ( 3 , I A ) ss SORTt 1 44E2 + . 25*TLAT ( 2 , I A )**2 ) 

22 CCNTINLE 

WRITE (LOUT. 34) ( I , ( TLAT < J , I ) . J= 1 ,3 ) , 1 = 1 .73 ) 

34 FORMAT ( / ( 5( l S,F5 . 0, 2F7 .2 ) ) ) 

DO 39 IA = I ,76 

7L0NG ( 2 , I A ) = 5*1 A - 1 00 

IF <TL0NG<2. IA) .GE . 180.) TL0NG(2,IA) = TLON G( 2 . 1 A >-360. 
TLDNG(l.IA) = TLONG (2«IA)+2.5 
29 CCNT INUE 

WRITE (L0LT.4I) ( I ,T LONG ( I * I ) « TLONG (2,1 ) ,1—1,7 6) 

41 FORMAT ( /( 6( I5.2F7.1 > ) ) 

DC 50 LAI = I ,73 

MLAT = 2-M0D(LA7,2) 

DO 49 LONG = 1 ,76 

JLATILCNG .LAT) - LAT 
* JLONG ( LONG , LAT) = LONG 

ALAT( LONG , LAT) = TLAT ( 1 , LAT ) 

ALONG (LONG. LAT) = TLONG( MLAT .LONG) 

49 CONTINUE 

50 CONTINUE 

WRITE (LOUT, 52) ( I . JLAT( I . 1) , JLONG ( I . 1 ) . ALAT < I ,1 ) , ALONG ( 1 , 1 ) 

$ .1=1,5548) 

2 FORMAT (/(5( 17.2 1 3 ,F 6 . I ,F7 • t ) ) ) 

WRITE (LL) ALAT 
WRITE (LU) ALONG 
RETURN 
END 
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SUBROUTINE BANDS 

cqmmo N 

*/ YPARAM/ LLAT.LLONG.LFNT1 » LPNT2 ,FREC( 1 A J .FREOI *FRE02 . SPACE < 2 ) 
% , WEDGE. DIRECC 20> ,DELT 

DATA LOUT/6/ 

NAMELIST 

S/L1ST02/ DIREC *FREO » CELT 
DC 61 IA * I .20 

D 1REC ( 1 A ) * IBAIA-S 
61 CONTINUE 

RAT « S.+* .OB 
FRECM 1 ) * FREOI 
DC 66 IA * It 13 

FREQUA+l) * RAT*FREO( IA) 

66 CCNT INUE 

*RIT£ (LOUT *L I ST 03 ) 

RETURN 

END 
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SLBRCLT I NE TRIG 
CCMMQN 

4/YPARAM/ LLAT »LLO NG , LFNT I .LPNT2 ,PREC{ 14 ) .FREOl .FRE02* SPACE (2 ) 

* , WEDGE. D1REC(20) .DEL T 
S/YSHORT/ TLAT(3.73) ,TL0NG<2.76) 

S/YYTRIG/ DLA 1 .COSS.D ANGLE »CLAT .TRAVEL, RAD.CTR.STR. 

S 0MIK,0I7.RAT(3.73) 

REALA E D L A T( 73 1 ,C0SS(20) . DANGLE (20 ) ,CLAT(2.73) .TRAVEL C 14 ) • 

* CTR (14)»STR(I4) .RA O.OMIN.OI7 
DATA LCUT/6/ 

REAL* 6 P 14 .SPEED 
PI4 s DATAN(IDO) 

017 * PI4/4SD0 
CM IN * P14/27D2 
RAD * 45D0/PJ4 

SPEED « 9 • 806D0/ 1 6D0 VFI 4 *3602/ 1 052DO 
C COMPETE TRAVEL (NAUTICAL MILES PER TIME STEP) 

DC 7E IA * I .14 

TRAVEL(IA) * SPEED *OELT/FREQ ( I A) 

CTR(IA) * DCOS( TRAVEL ( IA) AOMIN) 

STRCIA) s DSINCTRAVELC IA)AOMIN) 

70 CONTINUE 

WRITE ( LOUT .80) TRAVEL 
60 FCRMAT ( *0TRAVEL'/(7Die.l I) ) 

C COMPUTE AREA RATIO 

DC eS IA * 1.72 

RATCl.IA + l) s* TLAT (2.1 A + l ) /TLAT ( 2 . 1 A ) 

RAT ( 2 * I A ) « I. 

RAT ( 3. I A ) * TLAT ( 2 . 1 A) /TLAT (2,1 A+ 1 ) 
eS CCNTINLE 

RAT( I , I > * 0. 

R AT ( 2 » 73) a* I. 

R AT ( 2 • 73 ) * 0. 

C COMPUTE TABLE OF DIRECTI CNS 
DC 94 IA * 1 .20 

DANGLE! IA) =* DIREC(IA) 

IF <DIREC( IA) .GT. 100.) DA NGLE ( I A ) * DANGLE( IA )-36D 1 
DANGLE ( I A) s DANGLE ( I A) *01 7 
COSS(IA) as DC OS (DA NGLE ( IA) ) 

94 CCNTINLE 

C COMPUTE TABLE OF DOUBLE PRECISION COLATITUDES 
DC 99 IA as 1,73 

DL AT ( I A ) as ( 90 , ~TL AT ( 1 , I A) ) 401 7 
CL AT (1 » I A) as DCOS( DLAT( IA) ) 

CLA T ( 2 , 1 A ) = DSIN( DLAT(IA) ) 

99 CCNTINLE 

RETURN 
END 
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SUBROUTINE CHJTBND (I) 

CO MHO K 

5/YPARAM/ LLAT .LLONG, LPNT1 »LPNT2 »FREQ( 14)»FREQ1 .FREQ 2 , SPACE (2 ) 
9 , WEDGE. DIREC( 20) ,DELT 

S/YSHORT/ TLAT (3 *73) . TLONG (2.76 ) 

S/YYOUTE/ JTABLE. TABLE(4, 20.73) 

S/YYTRIG/ DLAT.COSS. DANGLE, CLAT. TRAVEL, RAD, CTR.STR. 

* 0MIN»017»RAT(3»73) 

I NTEGER42 JTABLE C 16 , 20 ,73) 

REAL*e DLAH73) .COSS (20) .DANGLE (20 ) .CLAT (2, 73) ,TRAVEL( 14). 

9 CTR( 14) o STR (14) ,RA D, ONI N, 017 
REAL* E SS * SMI 4 ) .DLAT2 »DAZ ,DTR .DLONG 
I N TEGER42 VERTEX ( 9 .4 ,E ) 

DATA LOUT/6/ 

DATA VERTEX/ 

9 1 . 1 .0 ,0 .1 .0.-76 ,-7 5.0 .0.1 .0 ,-l ,0 ,0,-t .-76.0. 1.0 . 1 . 2*- 1 , 

* 0.-77.- I .-76,0.1 ,0 2 ,2*— l ,-2.-7 7,- l, 

S l.l.C.l , 0.0. -75 .-76,0 ,0.1 .0. 1 ,1,0,1 ,-75.0. 1,0, 1.2. 1.1. 

* -74,1,-75,0.1,0.2,2.1.2,-74,1, 

* 2*— 1 .0,0,1 . 0,76 .77 ,0,0,-| .0 ,~1 ,0 .0,-1 ,76 ,0 ,- 1, 0 ,3*-l,0. 

* 75,- 1, 76,0.- 1. 0,-2 ,2*— 1 .-2. 75,-1 • 

S 2*- 1,0, 1,0, 0.77, 76 ,0,0, -1.0.1, 0,0,1 ,77.0 ,-1,0, -1.2, 1.1, 

S> 78, 1 ,77, 0,-1 ,0,2,2 ,1 ,2,78, 1 . 

$ 1,1,0,— 1,0,0, — 77,— 76, 0,0, 1,0, 2 *— 1 ,0,-1 ,—77 ,0,1, 0,1 , — 2, 

* 2*- l ,-78.-1 ,-7 7, C. 1, 0.24-2. -I , -2 .-78,-1, 

S 1 ,1, 0,0, — 1,0, — 76,— 77,0,0,1 ,0,1 , 0,0,1 , —76,0 ,1 ,0, 3*1,0, 

9 — 7£, 1,-76, C,1 ,0, 2,1,1 .2 .—75 , 1 , 

* 2*- 1 ,0,-1 .0 ,0,75,76 ,0.0 ,- t ,0 ,2*-l .0 ,-l ,75. 0,-1 , 0,- 1 ,-2, 

* 2*- 1, 74,-1 ,75. 0,-1 , 0 ,2*— 2 ,— 1 , -2. 74,-1 , 

S 2*-l ,0, 0,-1 ,0,76.75,0,0.-1 ,0,1 .0,0,1 ,76. 0,-1 ,0,-1, 1,1,0, 

9 77,1.76,0.-1 ,0.2,1 ,1 .2,77,1/ 

KFREC = I 

DTR = TRAVEL ( KFREC ) * CMI N 
RADIUS = TRAVEL( KFREC) 

WRITE (LOUT. 107) KFR EQ , FREQ ( KFRE C ) , R AD IUS 
107 FORMAT ( 1 H l • I 2 ,F 1 O. 7 ,F 1 0. 4 ) 

DO 110 IA = l ,5840 
TABLE( IA .1 . 1 ) = 0. 

110 CONTINUE 

DO 112 I A = 1 ,23360 
JTABLE ( I A , 1 ,1 ) -■ O • 

113 CONTINUE 

DO 2 OS LAT = 1 .73 

KWAD1 = 4* MOD (LAT, 2 ) 

WRITE (LOUT, 116) LAT ,TLAT( 1 ,LAT) 

116 FORMAT ( • OL AT =*,I3.* LATITUDE =*,F6.1/) 

DO 202 IDIREC =1,20 

KVAD2 = KWAD1 +1 D IREC/11 
C SOLVE SPHERICAL TRIANGLE 

DLAT2 =DARCOS (CTR (KFREQ) *CLAT ( 1 ,LAT) + 

* ST R ( KF REQ ) * CL A T ( 2 . LAT ) +COSS ( I DIREC)) 

S £ = .5D0*(DLAT( LAT) +DTR + DLAT2) 

SM( 1 ) = DSIN( SS) 

SMI 2) = D SI N ( SS— DLA T ( LAT ) ) 

S M ( 3 ) = DSIN(SS-DTR) 

SM ( 4 ) = DSI N( SS-0LAT2) 

DLONG = 2D0*D ATA N(DSQRT( SM (2)*SM(4)/(SM(1)*SM(3)) >) 

DAZ = 2D0*DATAN( DSQRT ( SM ( 3 ) *S M ( 4 ) / ( SM ( 1 )*SK(2 ) ) ) ) 

C INTERPOLATING FACTORS FOR LATITUDE 
FLAT = DLAT2— DL A T ( LAT ) 

IF (FLAT .G T • 0. ) KWAD2 = KWAD2+2 
FLAT = ABS(FLAT) *28.64789 
IF (FLAT ,GT. 1. ) STOF 130 
C INTERPOLATING FACTORS FOR LONGITUDE 
FLONG = DLONG 


IF (FLCNG .GT. A.) STOP 133 
FI a . 5*<FL0NG— FLAT ) 

F 2 a . 5* ! FLQNG+F LAT ) 

C INTERPOLATING FACTORS FOR AZIMUTH 
FA2 — DAZ 

IF (IOIREC .GT. 10) FAZ a -FAZ 
FAZ a 180.-57.29578*FAZ 
FAZ = FAZ/ 1 0. 

KAZ = FAZ 

GAZ = FAZ-FLOAT! KAZI-.5 

IF (GAZ .LT. 0.) GO TC IA9 

FAZ1 a l.-GAZ 

FAZ2 = GAZ 

KAZ 1 a 1 +KAZ 

KAZ2 a 2+ KAZ 

IF ( KAZ2 .EG. 21) KAZ2 a 1 
GO TO .155 
149 CONTINUE 

FAZ1 = l.+GAZ 
FAZ2 a -GAZ 
K A Z 1 a 1+KAZ 

IF^KAZZJ^EQ. 0) KAZ2 = 20 
155 CONTINUE 


IF 

(FI 

.LT. 

0. ) 

GO 

TO 

160 

IF 

(FI 

• GT. 

1.) 

GO 

TO 

175 

IF 

( F2 

.LT. 

1 • ) 

GO 

TO 

165 


GO TO 170 
C TRIANGLE I 
160 ITR = 1 

TABLE ( l.IDIREC.LAT) = F2 
TABLE (2, IDIREC.LAT) = -FI 
TABLE ( 3 » ID I REC »L AT ) = 1 .-FLAT 
HI TO 179 
C TR I ANGLE I I 
165 ITR = 2 

TABLE ( l . IDIREC , LAT) a FI 
TABLE! 2. ID1REC .LAT) a FLAT 
TABLE(3. IDIREC.LAT) = 1.-F2 
GO TO 179 
C TRIANCLE III 
170 ITR = 3 

TABLE! 1 . IDIREC .L AT) = F2-1 . 

TABLE! 2. IDIREC.LAT) a l ■ —FLAT 
TABLEI3. IDIREC.LAT) = 1.-F1 
GO TO 179 
C TRIANGLE IV 
175 ITR = A 

TABLE! 1, IDIREC.LAT) = Fl-1. 

TABLE!2. IDIREC.LAT) a FLAT 
TABLEI3, IDIREC.LAT) = 2.-F2 
179 DC 1 E3 1 A = 1.3 

JTABLE! 4*1 A-3, IDIREC.LAT) = 

S VERTEX! IA .1 TF .KWAD2+1 ) +LAT 

JTABLE (4*IA-2. IDIREC.LAT) = VERTEX ( I A+6 , ITR .KWAD2+1) 
JTABLE(4*IA. IDIREC.LAT) a VERTEX ('. A+3 , ITR. KWAD2+ 1 ) 
JTABLE(4*I A-l , IDIREC.LAT) a KAZI 
183 CONTINUE 

C TURN FRACTION OF ENERGY INTO BAND KAZ2 
DO 194 IA = 1 .3 

GAZ = TABLEdA .IDIREC, LAT)-FAZ2 
IF (GAZ .LT. 0.) GO TO 192 
TABLE! IA. IDIREC.LAT) = GAZ 
TABLE{4 .IDIREC .LAT) = FAZ2 
JTABLE! 13, IDIREC.LAT) = JT ABLE (4* I A-3 , I D I REC ,LAT ) 
JTABLE! 14 , IDIREC , LAT) a JTABLE ( 4* I A-2 . ID IREC ,L AT ) 
JTABLE! 16, IDIR EC .LAT) = JTABLE ( 4* I A, I D IREC . LAT ) 
JTABLE! 15, IDIREC. LAT) = KAZ2 
GO TO 195 

192 JTABLE!4*I A-l , IDIREC.LAT) = KAZ2 

PAZ2 = -GAZ 
194 CONTINUE 
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C ADJUST ENERGY FOR CONVERGENCE OF MERIDIANS 
DO 200 IA a I ,4 

KA s JTABLE(4* IA-3.IDI REC,LAT)-LAT 
IF (KA .NE. 0) TABLE! IA. ID1REC. LAT) = 

* TABLE ( IA , IDI REC.LAT)*RAT (KA+2.LAT) 

200 CONTINUE 

202 CONTINUE 

IF(LAT.E0.73> WRITE (LOUT. 204) ( ( JTABLE ( 4# I— 3 * J . LAT ) * 

S JTABLE (4* 1-2. J, LAT) . JTABLE (4* I- 1 » J ,LAT ) . 

* JTABLE (4* I ,J , LAT ) .TABLE ( I , J.LAT) . 1—1.4) • J— 1*20) 

204 FORK AT ( 4( 1 5 , I 4 .2 I 3 .F9 . 6) ) 

205 CONTINUE 
RETURN 
END 


DO 


23 5 


23 8 


241 


25 4 

255 

256 


26 2 
263 
26 4 


SL8RCLT I ME INBND (IN) 

CCMMON 

S/TYOLTB/ J TABLE , T ABLE (4.20.73) 

%// TAB2(6 ,20 .73) , MOOD, JTAB2 

INT6GER42 JTABLE ( 16.20 ,73 ) , JTAB2( 1 3 ,20,73) , M00D(3. 25) 
DATA LOLT/6/ 

LL « IN 
REA IND LL‘ 

DC 235 I A = 1 .25 

MOOD ( 1 . I A) = —2 1 28 C 
MOOD (2. I A) - 0 
MOOD ( 3 , 1 A) = 4-2 1280 
CCNT IN UK 

DO 236 I A = 1 ,0760 
TAB2( I A , 1,1) = 0. 

CCNT INUE 

DC 24 1 I A “ 1 ,10960 
JTAB2C I A • 1 ii l ) = 0 
CONTINUE 

DO 256 LAT = 1 ,73 

DO 255 IDIREC = 1,20 
DO 254 I A = 1.4 

KLAT — JTABLE ( 4*1 A-3. IDIREC. LAT) 

IF (KLAT -EQ. 0) GO TO 254 
IF (KLAT .EC. 74) GO TO 254 
KDIREC = JTABLE(4*I A-l . IDIREC , LAT) 

KOUNT = JTAB2( 1 ,KD IREC .KLATI+ 1 
IF (KOUNT .EQ. 7) STOP 250 
JTAB2(2*K0UNT. KDIREC. KLAT) = 

$ JTABLE(4*IA-2, IDIREC. LAT) 

JTAB2( 2*KOUNT+ 1 *KDI REC , KLAT ) = 

S MOOD ( LAT , 1 >+ 280*JTABLE(4*I A . IDIREC, LAT ) + 

* 1 4*( 1DIREC-KD IREC) 

TAB2(KCLNT .KOI REC, KLAT) = TABLE ( I A . 1 D IREC, L AT ) 
JTAB2( 1 .KDIREC .KLAT) = KOUNT 
CONTINUE 
CONT INUE 
CCNT INUE 

DC 263 LAT = 1 ,73 
ARITE (LU) 

$ ( JT AB 2 ( I . 1 .LAT ) . 1 = 1 ,260) , ( TAB2 ( I . 1 .LAT ) . 1=1, 1 20) 

IF( LAT .HE. 73) GO TO 263 
DO 262 IDIREC = 1.20 

KCLNT = JTAB2 ( l . ID IREC. LAT) 

ARITE (LOUT. 264) KOUNT, 

S ( JTAB2( 2*I . IDIREC.LAT) ,JTAB2( 2* 1+ 1 . I D I REC. L AT ) , 

* TAB2( I . IDIREC. LAT ) .1=1 .KOUNT) 

CONTINUE 

CCNT INUE 

FORMAT ( I2,6( 15, I7.F5.6)) 

REA IND LU 

RETURN 

END 


uuu 


S LBROLT I NE TR ANSP ( 1 OLD ,NEW) 

CC MMON/Y SHORT/TL AT ( 3 ,73) .TLONGI2.76) 

S//JTAB3! 13,2 0) .TAB31 6,20) ,JTAB4( 13, 14,20,73) ,T AB4 (6, 14, 20, 73) 
I NTEGER* 2 JTAB3 , JTAB 4 
DATA LOUT/6/ 

264 FORMAT (1 2,6(15, 17 »F9 ,6) ) 

LU=)OLD 

LL2sNEH 

DO 293 LA T = 1 ,73 

DC 2e5 1FREQ= 1,14 

LL 1— LL+IFREO 

READ ( LU1 ) JTAB3.TAB3 

DC 284 1 D I REC = 1 ,20 

DC 27e IA-1,13 

J 1 AB 4 ( I A , 1FREQ , I DIRE C »LAT )= JT AB3 ( I A , 1 D1REC ) 

278 CONTINUE 

DC 282 I A= 1,6 

TAB 4! I A, IFREO , ID 1 REC ,LAT) =TAB3 ( I A , I D1REC) 

283 CONTINUE 

284 CONTINUE 

285 CONTINUE 

WR1TEIL0UT ,294) L AT , TLAT ( 1 ,LAT ) 

DO 292 I A= 1 , 280 
KCLNT— JTAE411 ,IA,1 »L AT ) 

WRITE! LOUT ,264) 

SKOLNT, ( JTAB4( 2*1 , I A , 1 ,LAT) ,JTAB4 (2*1+1, 1A, 1 , LAT) , 

ST AB4 ( I , 1 A , 1 , LAT) • I = l .KCUNT) 

292 CONTINUE 

293 CONTINUE 

294 FORMAT!' TRANSPOSED PROPAGATE TABLE, LAT=',13, • LATITUDE =• , 
S F6.1/) 

WR I TE ( LU2 ) JTAB4 
WRI IE ( LU2 ) TAB4 
STOP 203 
END 


BIO 


PROGRAM ICEDECK 


I NTEGER*2 LANSE A ( 76.73) .LAND 1 ( 76 .35 ) .LAND2 ( 76 . 20 ) . LANDS ( 76 . 1 0 > 
ECU I VALENCE 

S (LANOH 1.1) .LANSEA (1.1) ) » ( LAND2( 1.1) .LANSEAC 1.36)) • 

£ (LAND3( 1 .1 ) .LANSEA (1 .56)) 

DATA LLAND.LOUT/l 1 .6/ 
lO FORMAT ( 2X , 212 »2X » 721 1.2X.2I1) 

DATA LAND 1/4 *2. 3*1 ,1 0*2.36*1 .23*2.2 3*2.27*1 .26*2. 

S 25*2,22* 1 .20*2.39*2, 1 ,1.2,1 .33*2 . 300*2 ,5*2 . I ,70*2, 

S 4*2, 1, 7 1*2, 4*2, 1 , l .70*2 ,4*2. 1 . 1 ,7 0*2 .4*2.1 . 1 .47*2, 1 ,22*2 , 

1, lc 70*2, 5*2,1 .1.69*2.4*2.3*1 ,69*2,5*2,3*1. 46*2,1,21*2. 
4*1,39*2, 1 . 1,27*2 .5*2.4*! ,39*2.3*1 ,25*2, 

1, 14*2.1, 1 , 17*2.3*1 ,2,4*1 ,26*2, 

5* 1, 13*2,3* 1 . 1 6*2 .0*1,26*2. 

4* l , 13*2,4* 1 , l 6*2 ,8*1 .26*2, 

5* 1, 12*2,4* 1, 1 6*2 ,8*1,26*2. 

5* 1 , 12*2.4* 1 , 1 6*2 .7*1.27*2, 

6* 1 , 13*2,3* 1 , 1 5*2 .7*1 . 27*2, 

6*1, 12*2.3*1,17*2,6*1,27*2, 

7*1, 10*2.5*1,1 7*2, 5*1 ,27*2. 

7* 1 , 10*2.6* 1 , 1 7*2,1 . 1 , 2. 1 ,20*2, 

8* l , 10*2,6* 1 , 1 7*2 , 1 , 1,29*2, 

7* 1 . 10*2.6* l . I e*2, 1,2,1 .28*2, 

8*1, 10*2.5*1,16*2,1 ,32*2, 

9* 1 , 10*2,5* 1 , 1 4*2 ,1 .5*2. I » 1 .27*2, 1 . 

9* 1 . 10*2,5* 1 ,20*2 .4*1 ,24*2, 1 , 

8* l, 10*2.6*1 , 1 2*2 . 1 ,3*2, 1 ,2 ,2 , 1 , 1 ,26*2, 1/ 


S 4*2 
£ 4*2 

*5*2,4* 
S 5*2 

5*2 
5*2 
5*2 
5*2 
5*2 
5*2 
4*2 
4*2 
4*2 
4*2 
3*2 
3*2 
3*2 


S 

S 

$ 

s 

s 

s 

s 

s 

* 

* 

* 

s 

$ 


DATA LAND 2/3*2. 8* 1,1 C*2 .7*1 . 14*2 . 1 , 1 , 2 , 1 , 29*2 , 1 , 




* 3 *2, 6* 1 , 12*2.7*1 .1 1*2 . 1 ,2 . 1 . 1 . 2 , 1 ,30*2 . 1 , 

S 4*2,5* 1 , 12*2,8* 1 , 1 0*2 . I ,2, 2, 1 ,33*2, 

£ 3*2,6* l , 12*2,8* 1 ,9*2 , 1 .1 ,2,2, 1 ,32*2, l , 

S 4*2,5*1.8*2,12*1,15*2,1,31*2, 

S 2*2.4* 1 ,9*2, 13* 1 , 1 4*2, 1,31*2,1. 

S 2.2, 1 ,2 ,3* 1 , 10*2, 12*1 . 1 0*2, 1 ,3*2, 1 ,29*2. 1,2, 

£ 1, 1 , 14*2, 1 2*1 .6*2, 1,5*2.1,33*2 ,1,1, 

S 2.1,1, 13*2. 13*1 .5*2, l ,1 ,3*2, 1 , 1 ,32*2,1 .1.2, 

£ 1, 15*2, 13*1 .5*2 . 1 , 4*2 . 1 . 1 , 30 *2 . 1 , 1 , 3* 2 , 

* 1.1, 14*2.15*1 ,3*2, 1. 1 .2,2, 3* 1 ,30*2,3*1 ,2 ,2 , 

£ 2. 1, 2, 1., 12*2, 1 0*1, 2,4*1 .3*2, 1 , l , 2 ,2 , 1 . 1 , 30 *2 , l , 1 .2 , 1 ,2, 1 , 

S 3*2, 1, 1 2*2. IS* 1 ,2,2 ,9*1 ,29*2. t .3*2.1 . 

S 3*2.1, 12*2, 10* 1 , 2, 3*1 ,2,2, 1 1 *1 ,27*2, 1 ,4*2. 1 , 

£ 17*2,9*1 .2.3* 1.2,1 2*1 .27*2,1 , l ,4*2. 

* 2,2,1, 13*2.27* l .28*2, l .2,2.1 .2. 

* 3*1.15*2.11*1,2,14*1.24*2.7*1.2, 

£ 3*1 , 14*2 ,7* 1 ,2. 2. 1 7*1 ,25*2,7*1 .2. 

S 4*1 , 14*2 .4*1 .5*2. 1 7*1 , 2 ,2, 1 , 2 1 *2 , 8* 1 , 

8 4* 1,15*2.1. 1.5*2, 3 *1.2. 13* 1, 2. 1,2, l ,2 0*2 ,9*1 / 

DATA LAND 3 /» 

S 4*1 , 14*2. 1 ,3*2. 1.2, 2. 1,1, 2, 1.2. 13*1. 2. 1.2, 2. 1,1 9*2 ,9*1 , 

£ 4*1 , 13*2. 1 , 1 .4*2,2 2*1 , 2 1*2, 1 0*1 , 

S 6*1 , 12*2 , 1 . 1 .2, 2.3*1 .2 .1 ,2 , 1 8*1 ,21*2 .9*1 , 

£ 5* 1, 14*2 .6*1 .2,2. 1 9*1 ,2,1 , 18*2. 10*1 . 

S 8*1,11*2,7*1.2, 20* 1.2,2,1,17*2 .9* 1 . 

S 6* 1 ,2. 1 . 10*2.30*1 , 18*2 , 1 0* 1*6*1 , 1 4*2.28* 1 , 18*2. 1 0* 1 , 

£ 8*1 . 10*2. 1,2,28*1, 2.2, 1 r 14*2. 11*1 , 

£ 8* 1.9* 2, 1,3*2.27*1 .3*2,1,13*2.11*1. 

S 1,1. 2, 4* 1,11 *2. 1,2, 2.1, 2.24* 1 ,3*2 , 1 , 1 ,12*2 , 1 0* 1 ,2, 1 , 

S 1. 1.2.3* l, 12*2. 1 ,3*2,1,2,24*1 .3*2,1 .12*2,10*1.2, 1, 

S 3*2, 1, 1 • 15*2,4* 1,2 ,25*1 ,2,1 ,6*2 ,14*1 ,3*2 , 1 , 

£ 1,2.2. 1 .5*2, 1 .1 . 9* 2, 3*1 ,2,31 *1 .3*2,15*1 ,2,2, 1, 

£ 1 , 2 . 2 , 1 . 1 , 4 * 2 , 1 . 1 , 3 * 2 ,1 , 1 , 5 * 2, 1 . 1 , 2, 1 , 1 , 2 , 28*1 . 4 * 2 , 14 * 1 . 2 , 2, 1 

£ 4 * 1 , 2 , 1 , 1 . 2 . 2 , 3 * 1 . 3 * 2 ,1 , 1 , 5 * 2 , 5*1 , 2 , 48*1 , 

£ 6* 1 .3*2 .5* l ,8*2.5* 1,2,2.25*1 ,4*2,18*1, 

£ 1 ,2, 4* 1 . 3*2. 6* l .9*2, 1 ,6*2, l . 2 , 1 8* 1 . 8*2. 8* 1 ,2 .3* 1 , 2 . 1 , 2 , l . 1 . 

£ 2. 4* 1.3 *2,7*1 . 14*2 ,1 ,3*2,1 2* 1 ,23*2,1 . I .3*2, 3*1/ 

WF ITE ( LOUT , l 0 -) LANSEA 


ORIGINAL Hmifc - 
OH POOR QUALITY 
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PROGRAM PAXXPAXX/MAIN 


COMMON 

•t/XXXXlO/ L W I ND * LL AND ,LPR OP , LOUT , LARCH 1 , LARCH2, LRS T 1 , LRST2 , 
$ LSCR 1 » LSCR2 • RUN I D 

S/XXTI ME/ K3TEPI , KSTE P2 , KSTEP3 , 1 ST EP . DELT , DELT AH, Y MPH , ZMDH, 
S lsk,ip.kot 

S/XSpGCl/ SPC)LU(2B0, 76,3) 
i>/ XSPEC2/ SPNEW( 280,76) 

$ / XP ARAM/ CA.CB.FREOl ,FRE02 , WEDGE, D I REC ( 20 ) , COS S( 2 , 1 0 ) , 

$ D E LOG, OMEGA (2 , 14) , FRE Q ( 14) , DOM ( 14). RBW ( 1 4 ) , 0 MM 4 ( 1 4 ) , 

S 0M3D0( 14 ) ,OMt 16( 14 ) , 0M25U4) .D0MM4( 14 ) , OMK (14) 

S/XXGRID/ L AT , LONG , KP C I NT 
S/XXLAND/ L AN SEA ( 7 6,7 3 ) 

S/XXWIND/ W l ND (2,76,73) 

J/XXXH 13/ FRUNT(4 ) ,Ht 3(76, 73) 

$/ XX PROP/ JPROP( 13 ,28 0 ,73) .PROP (6 ,2 8 0,73) 

$ / XXP ACK/ EEE(76) , SPE C 1 6 ( 2 80 , 76 ) , 

INTEGERS SPEC 16 

INTEGERS DELTAH.RUN ID. YMDH.ZMOH 
INTEGER+2 LANSEA.JPROP 
CALL WWQRK 

STOP ' • 

END 
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SLBRCLTINE WWORK 
COMMON 

S/XXXXIO/ LWIND. LLAND .LPROP, LOUT. LARCH1 .LAPCH2, LRST1 , LF3ST2, 

* LSCR 1 * LSCR 2 » RUNI D 

S/XXTIME/ KSTEP1 , KSTE F2 , KSTEP3 , I STEP ,DELT ♦ DELTAH. Y MDH , ZMDH, 

5 LSKIP, KBT 

S / XP ARAM/ CA »CB .FREG1 ,FREQ2 , WEDGE , D I REC( 20) , COS SI 2, 10), 

6 DELOG, OMEGA ( 2 ,14) . FREOl 14) ,DQM< 14) , RBW( I 4) ,0MM4 ( 14) , 

& 0M3DDI 14 ) ,OMt 16(14 ) ,0M25< 14) ,D0MM4( 14 ) ,OMK( 1 4) 

S/XXLAND/ LAN SEA (76,73) 

S/XXWIND/ WIND( 2,76,73) 

S/XXFROF/ JPROP( 13,280.73) ,PR0P(6 ,200,73) 

INTEGER* 4 DELTAH , RUN ID . YMDH , ZMDH i 

INTEGERS LANSEA, JPROP 

NAMELIST /WHAT/ L WI ND .LLAND , LPROP , LOUT , LARCH 1 , LARCH2 , 

S LRST 1.LRST2.LSCR1 , L5CR2 .KSTCP1 , KSTE P2 , KSTEP3 ,Y MDH. LSKIP . KDT 

Si CA , CB , RUN1D 
DATA LIN/5/ 

LXIND ~ 10 
LLAND = 11 
LPROP = 12 
LOUT = 6 
LARCH 1 a 13 
LARCH2 = 14 

LRST I a 15 
LRST 2 =16 
L SCR 1 = 17 
LSCR 2 =19 
KSTEF1 = 1 
KSTEF2 = 0 
KSTEF3 = 32767 
DELT = 5400. 

DELTAH = 3 
YMDH = ~ 1 
LSKIP = 0 
KBT = 1 
CA = 1.36E-9 
CB = . 3 066 
FREQ 1 = .04 
. 20 
i 8 • 

(LIN , WHAT) 

(LOUT, WHAT) 


FRE02 
WEDGE 
READ 
WR1 TE 
IF 


(LWIND .LE. 0) STOP 354 
IF (LLAND .LE. 0) STOP 355 
IF (LFRGF .LE. 0) STOP 356 
IP (LOUT .LE. 0) STOP 357 
IF (LOUT .EG. LW1ND) STOP 357 
IF (LOUT .EG. LLAND) STOP 357 
IF (LOUT .EQ. LPROP) STOP 357 
IF ( LARCH 1 .LE. 0) STOP 360 
IF (LARCH 1 .EQ. LWIND) STO° 360 
IF ( LARCH 1 .EQ. LLAND) STOP 360 
IF ( LARCH 1 .EG. LPROP) STOP 360 
IF ( L ARCH I .EQ. LOUT) STOP 360 
IF ( LARCH 2 .LE. 0) STOP 361 
IF ( LARCH 2 .EQ. LWIND) STOP 361 

if ( larch 2 .eo. lland) stop 36i 

IF ( LARCH2 .EO. LPROP) STOP 361 
IF (LAP.CH2 .EO. LOUT) STOP 361 
IF ( LRST 1 .LE.' 0) STOP 362 

IF ( LRST 1 .EO. LWIND) STOP 362 
IF ( LRST 1 .EQ. LLAND) STOP 362 
IF ( LRST 1 .EG., LPROP ) STOP 362 

IF ( LRST 1 .EQ. LOUT) STOP 362 
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IF (LhSTl .EO. L ARCH 1 ) STOP 362 
IF (LRST1 .EO. LARCH 2) STOP 362 
IF (LRST2 .LE. 0) STOP 363 
IF (LRST2 -EO. LWIND ) STOP 363 
IF CLRST2 .EO. ’LLANO ) STOP 363 
IF ILRST2 .EO. LPROP > STOP 363 
IF (LF.ST2 .EO. LOUT) STOP 363 
IF CLRST2 .EO. LARCH l ) STOP 363 
IF (LRST2 .EO. LARCH 2 ) STOP 363 
IF (LSCRl .LE. 0) STOP 36A 

IF CLSCR1 .EO. LWIND) STOP 36A 

IF (LSCRl .EO. LLANO) STOP 36A 

IF (LSCRl .EC. LPROP) STOP 364 

IF (LSCRl .EO. LOUT) STOP 364 j 

IF (LSCRl .EO. L ARCH 1 ) STOP 364 I 

IF ( L SCR ) .EO. LARCH2) STOP 364 j 

IF ( LSCR2 -LE. 0) STOP 365 
IF (LSCR2 .EO. LWIND) STOP 365 

IF ( LSCR2 .EO. LLAND) STOP 365 j 

IF (LSCR2 .EC. LPPOP) STOP 365 * 

IF ( LSCR2 -EO. LOUT) STOP 365 

IF (LSCR2 .EQ. LARCH 1 ) STOP 365 

IF ( L SCR 2 .EO. LARCH2) STOP 365 

IF ( L SCR 2 .EO. LRST1) STOP 365 

IF (LSCR2 .EQ. LRST2) STOP 365 

IF (LSCR2 .EO. LSCRl) STOP 365 

IF ( KSTEP) .LE. 0) STOP 366 

IF (KS7EP? .LT. KSTEP1) STOP 367 

IF (KSTEP2 .GT. 32767) STOP 367 

IF (KSTEP3 .LE. 0) STOP 370 j 

IF ( V MDH .LE. 0) STOP 371 
IF <MOD( YMDH, 1 00) .GT . 23) STOP 371 
KA = MOD ( YMDH .10000) 

IF (KA .EQ. 0) STOP 371 

IF (KA .GT. 3123) STOP 371 • 

KA = MODI YMDH , 1 OOOOO 0) 

IF (KA .EQ. 0) STOP 271 I 

IF < KA .GT. 123123) STOP 371 

IF (YMDH .GT. 99123123) STOP 371 

IF (LSKIP .LT. 0) STOP 372 

IF (LSKJP *GC. 32767) STOP 372 

IF ( KB T .LT. 0) STOP 373 

IF (KBT .GT. I) STOP 373 

IF ( C A .LE. 0) STOP 374 

IF (CE .LE. 0) STOP 375 

CALL BREW (LRST1) ! 

C ACQUIRE LAND-SEA TABLE [ 

CALL BREW (LLAND) 1 

CALL B RE AD 2 ( LLAND ,L ANSE A .5548 ) i 

C ACOUIFE INITIAL WIND FIELD • 

IF (LWIND .NE. LLAND) CALL BREW (LWIND) I 

WSKIP = LSKIP+1 ‘ 

DC- 381 I A = ' 1 . MSK I P • 1 

CALL BREAD4 ( LW I ND . W I ND , 1 1 096 ) 

381 CONTINUE 
C SET DATE AND TIME 
2 MDH = YMDH 

IF (KSTEP1 .EQ. 1) GO TO 387 
DC 386 I A = 2.KSTEP1 
CALL BUMP 

386 CONTINUE j 

387 CALL LODTAB ; 

C ACOUIFE PROPAGATE TABLES 

IF (LPROP .NE. LLAND) CALL BREW (LPROP) 

CALL B RE AD 2 ( LPR OP , J FROF . 265 7 2 0 ) 

CALL BREAD4 ( LPR OP . P ROP . 1 2264 0 ) 

DC 396. IA = KSTEPi .KSTEF2 
I ST CP = I A 
CALL WWAX1 
CALL BUMP 
CALL 0R*AX2 
396 CONTINUE 

RETURN 
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SUBROUTINE LODTAB 
COMMON 

. S/XXTIME/ KSTEP1 « KSTE P2 , KSTEP3 * I STEP,. DELT , DELTA!" ,Y MDH ? ZMDH, 
• * LSKIP.KBT 

S/XP ARAM/ CA.CB.FREQl ,I“R&02 , WEDGE , DI REC C 20 ) , COSS( 2 , I 0 ) , 

S DELOG, OMEGA (2,14) ,FREQ( 14) ,DOM( 1 4 ) , RBW ( 1 4 ) , OMM4 ( 14) . 

9 OM3DDI 14 ) , OM1 16(14) ?QM25( 14) ,DOMM4( 14 ) , OMK < 1 4) 

NAMELIST /TABS/D I REC ,C0S5 . OMEGA »FRE C, DOM »RBW , OMM4 , 0M3DD, 
S0M1 16,0M2S,D0MM4,0MK 
DO 4 15 I A - 1.20 

D I REC ( I A ) = 18*1 A- 9 
419 CONTINUE 

DC 425 IA = 1.10 

BB - .2141 592654FLCATI IA-1 ) 

COSSC 1 ,1 A) ~ COS(BB) 

C0SM2 , 1 A) ss SIN(BB) 

425 CONTINUE 

DELOG = .08*ALOG( 5. ) ' 

OMEGA! 1.1) “ .015625 > 

OMEGA (2,1) = 6.28218534FREQ1 
BE = EXP ( . 54DEL0G ) 

DO 422 1A = 3,28 

GMEGA ( I A , 1 ) = B B*0 MEG A ( 1 A- 1 , l ) 

432 CONTINUE 

DO 435 I A « 1,14 

FREQ ( 1 A ) * 0MEGAI2 , IA)/6. 283 1 853 
435 CONTINUE 

DO 436 IA - 2,13 

DOM(IA) = 0MEGA( 1 , IA + 1 )-OMEGA( l , I A) 

438 CONTINUE 

DOM ( 1 ) = DOM( 2) 

DCM(14) = . 25 /OMEGA ( 1 *1 4 > **4*OMEGA < 2 , »4)**5 
BACA « .3141 592654CA 
DC 449 l A = 1.14 

RBW(IA) = 6.283 1 653/DOM( I A) 

CMM4(TA) = OMEGAt 1 .1 A)**(-4 ) 

0M3DD ( I A ) = BACA+QM.EGA (2.1 A) **3*D0M< IA) 

□ M116( IA) = OMEGA ( 2 » I A ) **1 .16 
0M25(1A) = 0MEGA<2 .1 A> + *(-2.5) 

OMK(IA) = OMEGA (2 , I A )**2/9 .806 
449 CONTINUE 

DO 452 1 A ~ 1 , 13 

D 0 M M 4 ( I A ) = 0MM4( l A1-0MM4 ( IA + 1 ) 

452 CONTINUE 

D CMM 4 ( 14 ) = 0MM4 ( 14) 

DO 456 I A .= 1,14 

D0MM4(IA) = D0MM4( IA)*24. 039409 
456 CONTINUE 
C 24.03540S s= 9.806**2/4 
WRI 7E ( 6 , T ABS 1 
RETURN 
END 
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£ LORCLTI NE WWAX1 ' « 

COMMON /XXXXJ 0/ LWIND »LLAND , LPROP .LOUT . LARCMI , LARCH2 • LRST 1 ,_RS ( 
S LSCR 1 »LSCR2 »RUNI D /XXTIME/ KSTEPt .KSTEP2 * KST EP-3 » I STEP, DELT • I 
S DELTAH.YMDH.2MDHtLSKIP.K0T/XXFACK/ EEE(76 ) , SPEC 1 6 ( 280 , 76) < 

*/>SPEC2/ SPNEW! 280,76) /XXGP.ID/ LAT.LONG.KROINT /XXLAND/ ) 

* LANSEA (76,73) /XXtflND/ WJ ND ( 2 ,76 , 73 ) ! 

INTEGER72 SPEC! 6 , LAN SEA 
LOGICAL DEAD ST 
DEADS.T st ( I STEP -EG. I) 

CALL BREW (LSCR1) 

CALL GREW (LSCR2) 

IF (KOD! 1STEP.2) .EG. 0 .AND. I STEP .NE. KSTEP1 ) 

S CALL BREAD* ( LW I ND , WI ND , I 1 096 ) 

DO 820 LLAT * 1 ,73 ! 

LAT = LLAT 

IF (DEADST) GO TO 811 j 

CALL BREAD* ( LSCR 1 ,EEE ,76 ) 

CALL 3READ2 ( L SCR 1 .SFEC1 6 , 2 1 280 ) J 

CALL UNPAOQ ] 

GO TO 813 

61 1 DO 812 I A = 1 ,2 1 28 0 

SFNE W ( IA . 1 ) = 0. 

812 CONTINUE 

813 DO 81* LLONG = 3,74 

LCNG « LLONG 

IF (LANSEA! LLONG .LLAT) .EG. 2) CALL CMPE27 

814 CONTINUE 
CALL e ACQ 
EEEU) = EEE ( 73 ) 

EEE ( 2 ) = EEE ( 74 ) 

EEE ( 75 ) = EEE ( 3 ) 

EEE! 76) = EEE ( 4 ) 

DO 815 IAs l , 560 j 

SPEC 16! IA,1 ) = SPEC16(1A,73) 

SPEC 16< IA ,75) s S D EC 16! I A , 3 ) I 

815 CONTINUE ! 

CALL 9RITE4 ( L SCR2 ,EEE .76 ) 

CALL BRITE2 ILSCR2 .SPEC 16, 21280 ) 

820 CCNTINUE j 

RETURN ! 

END . * * 



B16 


SUBROUTINE CMPE27 

C REFERENCE SPECTRUM BASED ON U(19.5)**A 
COMMON 

$/> XT I ME/ K STEP l , KSTE R2 , KSTEP3 , I STEP .DELT .DELTAH, YMDH . ZMDH, 

* LSKIP.KBT 

S/XSPEC2/ SPNEWU 4*20 ,76) • 

S/XPARAM/ CA,CB,FRE01 ,FRE02 , WEDGE, D I REC ( 20 ) , COS S< 2 ,-l 0 ) , 

* DELOG, OMEGA ( 2. 1 4) ♦ FREO( 14) ,DOM( 1 4 ), RBW( 14) ,QMM4 (14 ) , 

S 0M3DD( 14J .0M1 16(14 ) ,OM2S(14) .DOMM4C14 ),0MK(1 4) 

S/XXGRIO/ LAT » LONG » KP C I NT 

a/XXLAND/ LAN SEA ( 76*72 ) 

S/XXWIND/ W 1ND (2, 76, 7 3) 

C S// DE, EX, OMSPEC. FMSPEC, SPREAD ,TPI G 
INTEGER*2 LANSEA 

REALA4 UXMUlT ( 66 ) ,UKAPPA<68) , GAMMAK (7S ) , DELT AK (75) 

DATA UXMULT/ 


S' 


& 

S 

X 

s 

s 

£ 

S 

s 

5 

* 

S 

s 


02245 51 30,. 032 SOS 74 9 
.033505060, .034030014 
• 0 3635 34 69, .036 5 66 276 
.03521 5548,. 03 9742 365 
. 04 1746046, .042224514 
. 04 4 0664 70,. 04 451 1437 
. 04 624 1 270 , . 046662 650 
.04 83! 3869,. 0487! 8 SI 9 
. 0S031 4064, .050707404 
.052263245, .052648445 
. 054 177510, .054557336 
. 056C69509, . 056446 177 
. 05754 9569. .05 8324 94 2 
.055826457, .060202 149 
DATA UKAPPA/ 


11959678, 

03457963, 


C2 199407, 1 


C l 01 9831 . 
.95722051 , 
.96374 1 52, 
.97012699, 
.95654205, 
.94305868, , 
.52570514, 

■ SI 648767, 

. S0240237, 
.85043847, 

. (7758294 , 
DATA GAMMAK/ 


. 08116600 
. 03 125 60,0 
.01576 14 0 
.00767146 
.99454 892 
,981 02220 
.0^740442 
, 95383593 
940377 17, 

. 92705081 
.91386 045., 
. 90080 024 
.88785915 
.87502354 


032705244 , 
034604220, 
037562479 , 
0402581 13, 
042695259, 
044951 047 . 
047080670 , 
049121290 , 
051 098861 
053032369 
054936363 
056822446, , 
, 058700256 , , 
060578078/ 


.032782464, .033070655, 
.035200157, .035801416, 
.038126069, .038676973, 
.040763505, .04 125 9258, 
.043158 832, .043615745, 
.045385666, *045015635, 
.047494955, ,047905950, 
.049521 158, .049918700. 


.05148857! , 
.053415122, 
.055314681 , 
.057198388, 
,059075572, 


.051076658, 

.053796804, 

.055692370, 

.057574074, 

.059450952, 


1.0596 1519, 1 
J . Or. 66564 2 . 1 
1 .01746383,1 
1 .00510420, 1 
.99186222, 

. 9782 9975 , 
.96468406, 
.951 1 3428, 


04698433, 

> n t 


01510145, 

00250289, 

98910370, 

97557563, 

96196653, 

94843744, 

9377 01 , .935 0302 2. . 9323 61 9 5, 
.92440198, .92175855, .91912054, 

.91123831, .90862127. 

.09820291, .89561024, 

.88528406, ,88271305, 

.87246776/ 


1.03938176, 
1, 0 24 ! 8266, 
1. 0 1267745, 
.99987330, 
.98645602, 
.972851 04 . 
.95925236, 
. 94574545, 


.90600933, 

.89302212, 

.88014607, 


£ . 1371 1 098, . 15359526 

S . 20213850, .20891603 
* . 2364284 3, .241 0401 6 

5 . 26125960, .2648532 2 

S . 281 17480, . 284 161 B 0 
S . 25600653, .3005053 1 
£ .31 269667, .31 497925 

S .32579964, .3278561 7 
£ .32767232, .33954982 

.34055516, .35029086 
. 2 Se 63 5"7 0, .36024603 
.36003 265, .3695401 2 
. 37 605 0 19, .3782 692 3 
DATA DELTAK/ • 


£ 

£ 

S 

s 

c 


. 75663491 , 
• 7 C2 23 0 35, 
. 68275549, 
. 67C01 786, 
. 66CB002 0, 


. 7406C86 1 
. 6988705 7 
.68027285 
. 66829403 
.65940474 


. 16647946 
.21517106 
.24543554 
. 26 832 56 6 
. 2870701 9 
.30310838 
.31722006 
.32908004 
.34 140100 
.35200083 
. 36183806 
.37103195 
.37967470/ 

. 72944406 
.69499353 
. 67795867 
.66665454 
.65821961 


.17721295, 
. 22099060, 
.24963734 , 
.27168624, 
.28990472 , 
. 30557857, 
.31 942090, 
.33187249. 


. 1 8640458, 
.22644076, 
.253 66 4 77, 


19472033, 
231 57303, 
2575341 7, 


.27494321, .27810388, 


.29266973, 
.30799847, 
.321 58342. 
.33383464, 


.34322672, .34502781, 

35368970. . 35535809, 

36341228. . 36496915, 
3725085 C, .373970 16, 


.2953691 4, 
.31037043, 
.3237091 8. 
.33576757, 
.34680504 , 
.35700656, 
.3665091 3, 
*37541727. 


.72090019. 
.69150484 , 
.67579202, 
.66509171 , 
.657001 17, 


.7082059 7, 
.68543228, 


.71399369, 

.68833522, 

.67375 55 8, « 671 63491 , 
.66359881 , .66217004, 
.65582620, .65469175, 
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5 .653S921B, .65253*12, .65150639, .65050999, .6*95*312, .6*860*13i 

St .6^7tPH8i • 6468037 0 » • 64 593972 » .6*50981 1. .6**2778*, .6*3*7780. 
S .6*269727, .6* 1 9351 2 ,. 6* 11 9060 , .6* 0*6292, .63975137, .63905526. 
£ .6363739*, .6377 068 1 ,. 6370533 2 ,. 636* l 293, .63578513. ,635169*6. 

$ .63*56 5*6, .6339727 2, .633390 81 , .63281938. .632 25806, • 63 1706*9. 

S .621 16*37, .63063139, .630107.13, .62959162, .62908*30, .62858501 . 
$ .626093*9, .62763952. .62713291 , .62666338, .62620077. .6257**09. 

S .62529554, .62*85252, .62**1 569/ 

REAl.** KAPPA, KAPPA1 

RFAL** EX( 20) ,DE ( 20) > 

REAL** TR1 G( 20,3 ) »QMSPEC( 1*) , PM5PEC (1**2), SPREAD ( 14,20,2) j 
ALPHA! X) * ■ 0081 /X** .23 

CrCS(X) a (3. 8663 239E-9* X- 1 . 5230871 E-* >*X* 1 . 

SI1N(X) * U .3496016E-1 1 *X-6 . 86096 1 6E-7 ) *X+ . 017*53293 
C N.B. THESE FORMULAS GOOD FOR ND .GE. 16 I 

C ACOUIRE (.STAR AND THETA SUB W 

X SPEED * WI ND < 1 s LONG ,LAT )— 2. 1 1 . * 

IF ( XSPEED .LT. 1.) GO TO 675 
KSPEED * XSPfeED 

C IF (KSPEED .EO. 0.) GO TO 675 

If (KSPEED .GE. 68) GO TO *9* 

XSPEEU * XSDEED-AINT (XSPEEO) 

Y SPEED a 1 »—X SPEED 
UX *5 X J NO ( 1 , LONG » l AT ) * 

$ ( YSFEED*UXMULT < KSPEED) +XSPEED*UXMULT ( KSPEED+ 1 ) ) 

KAPPA1 ss ( YSPEED*UKAPPA(KSPEED)+XSPEED*UKAPPA< K5PEED+1))/ 

* WINDfl , LCNG.LAT) 

GO TO *96 

*9* UX a XTND( 1 . LONG , LAT ) *UXMULT ( 68) 

KAPPA1 a UKAPPA ( 68) / XI ND ( 1 , LONG »LAT ) 1 

*96 TFX at XI ND< 2. LONG .LATJ + 1B0. | 

IF (THX .GE. 360.) THX a THV.-360. 

C WRITE (6, 10) CA.CB.FREQl ,FRE02, WEDGE . < DIREC ( I ) , I« 1 ,20) 

C 010 FCRMAT(/1X, *CA=* ,F10.2,2X, *CB=* .F10.2.2X, 'FREOla* .F10.2.2X, I 
C » *FREC2s» ,F 10. 2 ,2X» • WEDGE* • . F* . 2/ 1 X , «P I RECa • , 10 (FI 0.2.2X) ) 

KTH a THW/WEDGE ‘ 

IF (KTH -EO. 20) KTH * 0 
C COMPUTE COS( THETA BAR - THETA W) , ETC. 

C DC 1 11=1,10 

€ X R I TE ( C . 2 ) IT,COSS( 1 ,IT) ,COSS(> 5 1T) 

C 002 FORMAT ( / 1 X , UTs' , I 2 , 2X . • COSS ( 1 , l T ) = • , P8 . 2 , 2X , • COSS<2, fT)s* , 

C SFS.2I 
C 001 CONTINUE 

DC £03 |A a 1,20 
TRIG( I A , 1) = 0. 

503 CONTINUE 
KA = KTH 

ANG a THW— DIREC ( KTH+ 1 ) 

A NG2 a ANG** 2 ' ; 

CAA a CO CS ( ANG2 ) 

SA a SIIN< ANG2)*ANG 

C WR1TE(6,3) KTH. ANG.AKG2 ,CAA, SA j 

C 003 FCRMAT(/lX,'KTHa* ,F8,2.2X,'ANGa* ,F8.2,2X,» ANG2 = * ,F8.2,2X. 

C £ *CAA = .» .F8.2.2X, 'SA=* ,F8.2) 

DC 513 1A = 1.10 
KA a KA+1 

IF (KA .EQ. 21) KA a i • 

TRIG(KA.l) = COSS(l.lA) *CA A+COSS (2,1A)*SA 
513 CONTINUE 

DC 4 ITwl ,20 

C WR 1 TE ( 6, 7 ) 1 T ,TR IG(IT,1) , TR IG( IT, 2) ,TRIG(IT,3) 

C 007 FORMAT! /I X , • IT- • , 1 2 , 2X , • TR1 G ( I T , 1 ) = • . FI 0 .2 . 2X , *TRIG( IT, 2)=', 

C SF 10.2.2X, *TRIG( IT.3) = ' ,F10.2) 

004 CONTINUE 
KA = 10 

DC 521 I A a 1 ,20 
KA a KA+l 

IF ( KA .EO. 21) KA = l 

IF (TRIGtIA.l ).EQ. 0.) TRIG(IA.l) a -TRIG(KA.I) 

TRIG ( I A , 2) = TRIG! IA,1) **2 

IF (TRIG(lA.l) .GT. 0.) TRIG(IA,2) = 1 .-TRI G ( IA , 2) 

521 CONTINUE 

C COMPUTE CM, EGA ZERO AND OMEGA PEAK 


OMO * (. 92748721 *9. 0C6)/WINDU .LONG. LAT) 

C » 927 41721 * .74**. 25 

QMPE AK * .94574101*0*0 
C COMPUTE TOTAL. ENERGY 
ED * 0. 

DO 530 I A * I »aao 

ED * ED + 5PNEW ( I A * 1 .LONG) 

530 CONTINUE 

UX4 * UX**4 
0*04 * OMO**4 

GUX4 * OM04/( .0001*. 25*9.806*9.806) 

UG * UX/9.806 

C COMPUTE REFERENCf SPECTRUM 
DC 540 JA a 2,14 
PMSPEC ( 1 A • 1 ) * 0. 

BA * QM04*0MM4 ( 1 A) 

IF (BA *LT . 00.) PKSPECUA.l) * EXP(-BA) 

54 0 CONTINUE 

BA » ( .25*9. 006*9.006) /QM04 
PMSPEC ( 1*2) * BA*PMSPEC( 2,1 ) 

D0S45IA®3»)4 

PMSPEC ( I A- 1 ,2 ) * 0 A* ( PMSPEC ( I A * l I —PMSPEC ( I A— 1,1)) 

54 5 CONTINUE 

PMSPEC (14,2) * BA*(1 .-PMSPEC ( 1 4 , 1 ) ) 

0MP25 * CMP£AK**2 • 5 
DC 572 IA a 1,14 
S = 0MP25*0M25( 1A) 

IP ( $ .GT. 1 . ) S a 1 ./S**2 
S « 1 5 » 0 04 96*S 
KS « S 

IP (KS • NE • 0) GO TO 559 
DO 657 IB = 1.20 
• TRIG ( IB, 3) » 0. 

IF (TRIG( IB. 1 ) .GT. 0.) TRIG(IB,3) a TR I G (, J B. 1 ) 

557 CONTINUE 

GO TO 555 

559 XS a S-AJWT(S) 

DC 6C4 IB » 1.20 
TRIG (1B*3) a Q, 

IF ( TP 1 G ( IB, 1 ) .GT. 0. ) 

$ TRIGUB.3) = ( 1. +T»IG< IB, 1) )**K5*( 1 .+XS*TRI6( IB, 1 ))- I. 

564 CONTINUE 

565 BA = Or 

DO 568 IB a 1,20 

BA a 13A + TRI 3 ( IB, 2) 

56 8 CONTINUE 

DO 571 IB a 1 ,20 

SFREADU A.IB.2) a PMSPEC ( I A , 2 ) *TRI G ( I B . 3 )/BA 

571 CONTINUE 

572 CGNUNLE 

C CYCLE THROUGH FREQUENCIES 
DC 674 I A a 1,14 
C CALCULATE ALPHA 
AL = 1 . 

IF(ED.GT.O.) ALaALPHA (ED*GUX4 ) 

OMSPECUA) a 0. 

DO 560 IB a |,20 

OMSPECUA) a OMS PEC ( 1 A ) 4-SPNEW ( I A , 1 B , LONG ) 

580 CONTINUE 

C WRI TE (6,20) AL.ED ,OM£PEC, (SPNEWU A, IB .LONG) , IB = 1, 20) 

C 020 F0FMAT(/1X, *ALa* ,F10 .1 ,2X , »EDa* ,F10 .1 ,2X, 

C * ,/lX, ’OMSFECa* .FI 0.1 ,/lX, ‘SPNEWa* , I 0< F10 . 1 ,2X ) ) 

ED a ED—OMSPEC ( 1 A) 

C CUT BACK HIGH FREQUENCIES 

CL T a AL*DQMM4 { I A) . ' 

IF (OMSPECUA) .LE. CUT) GO TO 589 
CLT = CUT/QMSP£C< I A) 

C WR J TE (6,21 ) ED. CUT *D CMM4 

C 02 1 FCRMAT(/1X, *ED=' , F10 .1 , 2X , 'CUT=' ,FIO. 1 ,2X, • D0MM4a>» ,FIO. 1 ) 

00 ses IB a 1,20 

SPNFW< IA , IB, LONG ) ‘a SPNE W ( IA, IB, LONG) *CUT 
588 CONTINUE 

C WRITE (6,22) CUT , ( SPNEW ( I A . I B . LONG ) , IBa J , 20 ) 


C 022 FORMAT!/ 1 X, *CUT=‘ ,F1 0. 1 ,2X. • SPNEWS* , 10 (F10 . 1 , IX) ) 

589 CONTINUE 

C DISSIPATE OPPOSING BANDS 

DO 599 I B = I » 1 0 » 

IF ( SPNE W ( I A . I B . LONG ) .EQ • 0.) GO TO 598 - 

IF ( SPNEW< I A , IB4-I0 .LONG) .EQ* 0.) GO TO 598 

BA - A BSC SPNE W( I A • I B » LONG) — SPNEtt ( I A » I B+l 0 « LONG ) )/ 

5 C SPNEWC I A .1 B .LONG) 4- SPNE W (IA*IB+10. LONG ) ) 

SFNEWC IA , IB, LONG ) = SPNEW! I A, I B,LCNG)*BA 
SFNEWC 1 A ,IB+1 O.LCNG) « SPNEW ( I A » I B+10 * LONG )*B A 
698 CONTINUE | 

C COMPUTE DOWNWIND ENERGY IN FREOUENCY BAND 
BA - 0. 

DO 606 IB - 1 ,20 

IF ( TRIG! IB * 1 ) ) 606,603.605 , „ , 

603 BA = BA+.5*SPNEW (I A , IB, LONG) 

GO TO 606 

605 BA = BA + SPNEWCIA , IB, LONG) 

606 CONTINUE 

C IF UNDERDEVELOPED GRO I* . IF OVERDEVELOPED, REDISTRIBUTE. 

IF (BA ,GE. AL*PMS PECCIA.2)) GO TO 647 

C GROV 

PSI = OMEGA (2»IA)*UG 
KAPPA = KAPPA1+0M1 16 (I A) 

XKAFPA - KAPPA* 1 00 .-I . ■ 

KKAPPA = XKAPPA 

IF (KKAPPA .GT. 0> GO TO 618 

GAMMA - .0027422197 

DELTA = .0151326985 

GO TO 622 

618 X KAPPA = XKAPPA-AI NT(XKAPPA) 

YKAFPA = 1. -XKAFPA 

GAMMA = KAPPA* ( YKA FPA*GAMMAK( KKAPPA ) +X KAPPA* GAMMAK (KKAPPA+ I 
DELTA = KAPPA* ( YKA FPA*DELTAK< KK APFA) +XKAPP A* DELT AK ( KKAPRA+ l 
622 ANUM = UX4*0M3DD( I A) 

DD 644 IB = 1,20 

IF ( TRIG ( IS , 1 ) .L E. 0.) GO TO 644 

. SPREAD! IA , IB . 1 ) = AL* SPREAD (I A, IB, 2) I 

IF ( SPNEWI IA ,1B, LONG) ,GE. SPRE AD ( I A , I B , 1 ) ) GO TO 644 j 

PSI1 = PSI *T RIG( IB,1) 

BX = CB*(PSI 1**2-. 0004) 

1FC.KBT.NE . 0. ) 

I EX = BX+CB/ <503 .3+ < 2042E3+1 2203E4*PS 1 1 )* (PSI I— .031 )**2 ) 

IF (BX » LE . 0.) GO TO 644 

ET = BX70MEGA (2, IA.) 

BDT = BT*DELT ' 

IF (BDT .LT. 77.) GO TO 637 i 

SFNEWC IA, IB, LONG ) = SFRE AD ( I A . I B . 1 ) ! 

GO TO 644 ! 

637 AT = ANUM/ ' . i 

S < ( GAMMA+( OMK< I A)*TRIG( IB, l )— K APPA ) ** 2/ GA MM A ) * • ! 

S (DELTA+CMK(IA)**2*TF. IG(IB»2) / DELTA ) ) 

AB = AT/BT 

• SPNEWdA.IB.LONG) = ( SPNEW ( I A . I B. LONG ) +AB) *EXP (BDT >- AB 
IF < SFNEW(I A , IB, LONG) .GT. SPRE AD ( I A , IB , 1 ) ) 

£ SPNE W( 1 A ,1 B .LONG) = SP RE AD ( I A , I B , 1 ) 

644 CONTINUE 

GO TO 669 

C REDISTRIBUTE OVER ANGLES 

647 SX = 0. • i 

£Y = 0 . * 

DO 660 I B = 1 , 20 ! 

E X ( I B ) = 0 . j 

DE( IB) = 0. ■ i 

IF ( TR IG( IB,1 ) • LE. 0.) GO TO 660 j 

EB = SPNEW( I A • I B i LONG ) — AL* SP RE AD ( IA, IB, 2) , 

IF (BB> 655,660,658 ! 
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655 DE ( I 5 ) = -B8 

SY = SY-BB 
GO TO 660 

658. EX(IB) = BB 

SX = SX+BB 
CONTINUE 


66 0 
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IF ( SX .£0. 0. ) GD TO 669 
IF ( SY .£0. 0.) GO TO 669 
RATI « SX/BA 
RAT 2 = RAT 1 ♦ SX/ SY 
DO 668 16 = I .20 

3 F ( Tr? 1 G ( I B » 1 ) . G T . 0.) 

S6 SPNEM I A, IB .LONG) = 

S S=>NE W( I A. I B »LCNG)~RAT1*EX( IB)+RAT2*DE< IB) 

668 CONTINUE 

669 GMSPEC ( I A) = 0. 

DO 672 IB = 1.20 

OMSPEC(IA) = CMS PEC < I A)+SPNEW ( I A. IB. LONG) 
672 CONTINUE 

ED = ED+CMSPEC ( I A ) 

67 4 CONTINUE 

675 CONTINUE 
C COMPUTE EPS 

ED = 0. 

DC 68C l A = 1 »?80 

ED = ED+SPNEW(I A.l .LONG) 

68 0 CONTINUE 

EFS = ED*GUX4 

RETURN 

END 
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5 LBR CUT X NE OR WAX 2 

COMMON / XXPACK/ EEE ( 76 ) , SPEC 1 6 ( 2 00 , 76 > / XXXXIO/ LW1ND.LLAND, LPROP . 

S LOUT * LARCH 1 • LARCH2 »LRST l .LRST2 »LSCR1 » LSCR2 « RUN I D /XXTIME/ KSTEPl . 
S , K STEP2 » KSTEP3 * I ST EP * DELT * DELT AH » YMDH »ZMDH ,LSK I P » KQT /XSPEC1/ , 
S SPOLDJ 280, 76.3) /XSPEC2/ SPNEW (280 , 76 ) / XXGR ID/ LAT , LONG , KPO INT 
S/XXLAND/ LANSEA< 76,72) /XXWI ND/ WINDi 2,76,73 ) /XXXHt 3/ JFR0NTC4), I 
S HI3(76.73) ./XXPROP/ J PROP ( I 3 , 200 , 73 ) , PROP (6 , 280 , 73 ) 

INTEGER** DELTAH , RUN ID • YMDH , ZMDH 
JNTEGER+2 LANSEA . JPRCP , MOOD ( 73 ) , SPEC16 

DATA MOOD/3.1 .2. 3, 1,2, 3. 1,2, 3, 1,2, 3. 1.2. 3. 1.2, 2, 1,2, 3, l, 2. 3, 1,2, 3. 

S I, 2, 3, 1,2, 2, 1,2. 3, 1,2, 3, 1,2, 3,1 ,2.3,1 ,2, 3, 1,2, 3, 1*2, 3, 1.2, 3, 1,2, 
£>3,1, 2, 3,1 ,2,3,1 ,2,3, 1,0, 0/ 

JFRONT ( 1 ) = RUNID 

JFRONT ( 2 ) = 1 STEP 

JFRONT (3 ) = 7 MDH 

JFRONT (A) = MOD( I STE P » KSTEP3 ) 

IF ( JFRONT (4 ) .EO. 0) CALL BRJTE4 ( L ARCH1 , JFRONT, 3 ) 

CALL BREW (LSCR1) 

CALL BREW (LSCR2) 

CALL BREAD4 ( LWI ND , W I ND ,1 1 096) 

DO 622 I A = 1 ,5548 
H 13 ( IA , 1 ) = 0. 

622 CONTINUE 

CALL GRAB ( 1) 

CALL GRAB < 2) 

DC 830 LLAT = 1,73 
LAT =; LLAT 
ML AT = MOOD (LLAT) 

CALL GRAB (MLAT) 

DO 623 I A = 1 ,21280 
SFNEWCIA.l) = 0. 

823 CONTINUE 
DO 627 LLONG = 3,74 

IF < LANSEAfLLONG ,LLAT) .NE. 2) GOTO 827 
LONG - LLONG 
KPC1NT = 76*LLAT +LLC NG 
DO 025 I A = 1 ,280 

KPART = JPROP( 1 , 1 A, LLAT) 

IF (KPART .EQ. 0) GO TO 825 
DO 824 IB = 1 . KPART 

KB = KPO I NT+ JPROP (2* I B t I A , LLAT ) 

IF ( LANSEA(KB~76, 1 ) . NE. 2) GO TO 024 
KC = IA+JPR0F( 2 + 1 B+l .1 A , LLAT ) +2 1 280 
SPNEWt T A .LLONG) - SPNE W < I A , LLONG ) + 

S. SPOLD ( KC .LLONG , 1 ) + PROP ( I B ,1 A. LLAT) 

824 CONTINUE 

825 CONTINUE 
CALL C M D E27 
EE = 0. 

DO 326 I A = 1 .28 0 

EE = EE+SPN6W( 1A .LLONG) 

826 CONTINUE 

' H12< LLONG .LLAT) = 4.*SQRT(EE) 

827 CONTINUE 
CALL PAOQ 

IF ( JFR0NT<4) .NE. 0) GO TO 028 

CALL BRITE4 ( LARCH 1 , EEE ,76 ) 

CALL BR1TE2 < LARCH 1 , SPEC1 6 ,2 1 200 ) 

828 CALL BRITE4 ( L SCR 1 ,EEE » 76 ) 

CALL 3RITE2 ( LSCR l , SPEC 1 6 , 2 1 280 ) 

83 0 CONTINUE 

CALL BRITE4 ( LARCH 2 , JFRONT , 5552) 

IF ( JFRONT (4 ) .EQ . 0) CALL BEND ( LARCH 1 ) 

WRITE ( 6, 82 1 ) JFRONT, (<I,J,H13(J,1), J=3 ,74) .1=1 ,73 ) 

831 FCRMAT(lHl,4I10//(9( 14, I3.F7.3) )) 
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SUBROUTINE GRAB (I) 

COMMON / XSPEC 1/ SPOLD < 28 0 ,76 ,3) /XSFEC2/ SPNEW <280 ,76) 
S/XXPACK/ EEE (76) , SPEC1 6(280, 76 ) 

COMMON 

S/XXXXIO/ LWIND.LLAND ,LPROP,LOUT,LARCHl ,LARCH2, LRST1 .LRST2. 

S LSCF 1 , LSCR2 , RUN I D 
I KTEGER*2 SPEC 16 
IF (1 .EQ. 0) RETURN 
CALL BREAD4 ( LSCR2 ,E EE ,76 ) 

CALL BREAD2 ( LSCR 2,S PEC 16 ,2 1 2 80 ) 

CALL UNPAQO 

DC 821 IA = 1,21280 

SPOLD ( I A , 1 , 1 ) = SP NEW ( I A , 1 ) * 

821 CONTINUE 
RETURN 
END 



80 0 

80 1 
802 

803 

805 
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SLBRCLTINE PAQQ 

COMMON / XSPEC2/ SFNE W <280 .76) /XXPACK/ EEE ( 280 ), S PEC I 6 ( 280 , 76 ) 
INTEGER*2 SPEC! 6 
DP 80S I A = 1 ,76 
EEE < IA ) = 0. 

DO 800 I B = 1 .280 

IF (SPNEWi IB. IA) .GT. EEE ( I A) ) EEE(IA) = SPNEWCIB.IA) 
CONTINUE 

IF (EEE(IA).GT. 0.) GO TO 802 
• DO 801 IB = 1 .2 80 • 

SPEC 16(IB»IA) = 0 
CONTINUE 

GO TO 805 % 

BA = 32767. /EEE l IA) 

DO 803 IB = l .280 

SPEC 16( IB.IA) ~ EA*SPNEW(IB.IA)+.5 
CONTINUE 
CONTINUE 
RETURN 
END 



806 

* 

eo7 

808 
81 0 


SLBRCLTI NE UNPAQQ . } 

COMMON / XSPEC2/ SFNE W < 2 80 ,76 ) /XXPACK/ EEE (280 ') • SPEC 1 6( 280 , 76 ) 
INTEGER*2 SPEC 1 6 , ' f 

DG 81 C I A = I ,76 1 

IF (EEE(IA).EO. 0.) GO TO 807 ! 

BB = EEE ( I A ) /3 2 767 « ! 

DO EC6 IB = I ,280 | 

BC=SPEC I 6 ( I B , I A ) 

SPNEW(IB.IA) = BB*BC ! 

CONTINUE J 

GO 1C 810 

do eca ib - i ,280 ; 

SPN5 ¥ ( IB ,IA) = 0 . 

CONTINUE 

CONTINUE I 

RETURN i 

END : 
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SUBROUTINE BREAD2 (LU, ARRAY, NPOINT) 
INTEGER*2 ARR AY( NPOI NT) 

READ ( LU ) ARRAY 

RETURN 

END 


B26 


StOROtT I N6 8PEADA <LU, ARRAY , NpDINT> 
INTEGER*^ ARRAYtNPOl NT) 

READ (UU) ARRAY 

RETURN 

END 


B27 


SLBRCLT I NE BR!TE2 <L U, A PR AY . NPOI NT ) 
INTEGER*2 ARR AY< NPOI NT) 

*RITE ( LU ) ARRAY 
RETURN „ 

END 
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SUBROUTINE BRITE4 (Lt, ARRAY, NPQINT) 
INTEGER** ARRAY! NPOI NT) 

HRITE (LU) ARRAY 

RETURN 

END 
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SUBROUTINE BREW ( UU) 

REWIND LU 

RETURN. 

/ END 



SUBROUTINE BEND < LU) 

ENDF1LE LU 

RETURN 

END 


APPENDIX C 
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